莫曉梅，涂新軍，2，3*，王 天，謝育廷，趙國(guó)羊

(1.中山大學(xué)土木工程學(xué)院水資源與環(huán)境研究中心，廣東 廣州 510275；2.廣東省華南地區(qū)水安全調(diào)控工程技術(shù)研究中心，廣東 廣州 510275；3.南方海洋科學(xué)與工程廣東省實(shí)驗(yàn)室，廣東 珠海 519000)

干旱是指由于水分的收支或供求不平衡形成的水分短缺現(xiàn)象[1]。 干旱在世界范圍內(nèi)仍是發(fā)生頻率最高、影響范圍最廣、致災(zāi)損失最大的自然災(zāi)害之一[2-3]。 全球因干旱導(dǎo)致的經(jīng)濟(jì)損失高達(dá)60 ～80億美金，遠(yuǎn)超其他自然災(zāi)害[4-5]。 無(wú)論是基于水循環(huán)組分角度或研究對(duì)象要素角度，均可以將干旱分為氣象干旱、水文干旱、農(nóng)業(yè)干旱和社會(huì)經(jīng)濟(jì)干旱4類[6]。 其中，氣象干旱主要是基于降水量不足而導(dǎo)致的干旱，水文干旱則是由于徑流量不足而導(dǎo)致的干旱[7]。 氣象干旱作為所有干旱發(fā)生的源頭，也是最普遍和最基本的干旱類型，水文干旱則更能反應(yīng)旱情，且影響因素較氣象干旱復(fù)雜，影響程度更深[8]。 二者反應(yīng)了干旱事件的不同階段，通常認(rèn)為水文干旱是氣象干旱的延續(xù)，并直接影響農(nóng)業(yè)干旱和社會(huì)經(jīng)濟(jì)干旱。 此外，從氣象干旱到水文干旱的過(guò)程會(huì)受到土地利用因素、產(chǎn)匯流過(guò)程以及土壤性質(zhì)等因素的影響[9-11]。

總體而言，研究基于不同水文氣象要素形成的氣象干旱到水文干旱的傳遞，可以在機(jī)理上揭示干旱傳播途徑和干旱形成過(guò)程。 因此，在干旱傳遞的過(guò)程中，量化傳遞規(guī)律是十分必要的。 有研究從干旱指數(shù)相關(guān)性的角度，分析水文干旱滯后于氣象干旱的時(shí)間，以表征干旱傳遞特征[12]。 除通過(guò)干旱指數(shù)判斷外，直接通過(guò)降雨、地下水和徑流等影響要素對(duì)干旱傳遞進(jìn)行判斷也是一種適合的判斷方法[13]。盡管上述2 種方法在一定程度上解決了定量化閾值的問(wèn)題，但都將2 種類型的干旱進(jìn)行了剝離，這將會(huì)導(dǎo)致2 種干旱無(wú)法在時(shí)間上統(tǒng)一。 同時(shí)，2 種不同干旱指數(shù)之間的分布狀況不一致，也會(huì)造成結(jié)果的失真。

因此，本文基于之前學(xué)者研究存在的問(wèn)題，在東江流域提出一種更適合定量化判斷干旱傳遞閾值的方法。 首先，計(jì)算東江流域月尺度的標(biāo)準(zhǔn)化降水指數(shù)(Standardized Precipitation Index，SPI)和標(biāo)準(zhǔn)化徑流指數(shù)(Standardized Runoff Index，SRI)指數(shù)。 其次，通過(guò)游程理論提取干旱事件及其特征，如歷時(shí)和烈度，并對(duì)發(fā)生時(shí)間部分或全部重疊的氣象干旱事件和水文干旱事件進(jìn)行匹配。 最后，利用二變量Copula函數(shù)，構(gòu)建氣象干旱和水文干旱的特征屬性的聯(lián)合分布模型，通過(guò)條件概率分析，識(shí)別氣象干旱向水文干旱傳遞閾值。

1 研究區(qū)域概況及數(shù)據(jù)來(lái)源

東江為珠江流域三大水系之一，流域控制站博羅站以上的集水面積為25 325 km2。 東江流域?qū)儆趤啛釒Ъ撅L(fēng)氣候，流域多年平均降水量為1 500 ～2 400 mm，降雨時(shí)間上分配不均，4—9 月為雨季，降水量達(dá)全年降水量的80%，10 月到翌年3 月為旱季，降水量?jī)H為全年降水量的20%。 東江擔(dān)負(fù)著為珠三角地區(qū)4 000 多萬(wàn)人口提供生產(chǎn)、生活、生態(tài)用水的重任，因此研究東江流域干旱情況對(duì)社會(huì)和經(jīng)濟(jì)發(fā)展有著重要意義。 研究數(shù)據(jù)為中國(guó)國(guó)家氣象局氣象信息中心提供的東江流域內(nèi)34 個(gè)主要?dú)庀笳?956 年4 月至2009 年3 月的月降水?dāng)?shù)據(jù)和博羅水文站1956 年4 月至2009 年3 月的實(shí)測(cè)月徑流數(shù)據(jù)，并根據(jù)流域三大水庫(kù)(新豐江水庫(kù)、楓樹(shù)壩水庫(kù)和白盆珠水庫(kù))的月徑流調(diào)蓄量，還原得到天然徑流數(shù)據(jù)。 站點(diǎn)分布見(jiàn)圖1。

圖1 東江流域及站點(diǎn)分布

2 研究方法

2.1 干旱指數(shù)計(jì)算

由于SPI指數(shù)可以表征一定時(shí)期內(nèi)降雨發(fā)生概率，該指標(biāo)不但適用于多時(shí)間尺度的氣象干旱描述，而且能夠很好地反映干旱的強(qiáng)度和持續(xù)時(shí)間，因此被廣泛使用[14]。 SRI指數(shù)與SPI指數(shù)類似，計(jì)算方法簡(jiǎn)單，同樣適用于多時(shí)間尺度分析，對(duì)資料缺乏、地形復(fù)雜的區(qū)域也有很好的適應(yīng)性，目前在水文干旱識(shí)別研究中應(yīng)用廣泛[15]。 因此，本文分別采用SPI指數(shù)和SRI指數(shù)來(lái)表征氣象干旱和水文干旱，2種指數(shù)的計(jì)算方法如下。

SPI指數(shù)計(jì)算采用常用3 參數(shù)Log-logistic概率分布模型，其累積概率函數(shù)見(jiàn)式(1):

式中 (x|Xi=1，…，n)——降水系列；n——樣本系列的長(zhǎng)度；α、β、γ——基于樣本系列通過(guò)線性矩法擬合而得的參數(shù)。

對(duì)累積概率進(jìn)行標(biāo)準(zhǔn)化轉(zhuǎn)換，有:

因此，干旱指數(shù)SPI為:

式中，C0=2.515 517；C1=0.802 853；C2=0.010 328；d1=1.432 788；d2=0.189 269；d3=0.001 308[16]。

以SPI指數(shù)的計(jì)算理論為基礎(chǔ)，以徑流值代替原SPI計(jì)算過(guò)程的降水值，即可得到標(biāo)準(zhǔn)化徑流指數(shù)SRI[17]。

2.2 干旱事件識(shí)別

游程理論是分析時(shí)間序列的一種方法，經(jīng)過(guò)長(zhǎng)時(shí)間的發(fā)展和應(yīng)用，目前被廣泛應(yīng)用于干旱事件的識(shí)別中[18-19]。 游程理論可以有效地從連續(xù)的干旱事件中提取干旱事件的歷時(shí)D(Duration)和烈度S(Severity)。 歷時(shí)是指干旱事件持續(xù)時(shí)間，烈度是干旱事件中指標(biāo)值與截取水平差值之和，研究中通常取其相反數(shù)。

本文不僅通過(guò)游程理論對(duì)干旱事件進(jìn)行識(shí)別，同時(shí)為了避免次要干旱事件導(dǎo)致干旱傳遞模型被破壞，同時(shí)采用了干旱事件的融合和剔除方法[20]。 具體步驟如下[21]:①當(dāng)指數(shù)值小于R1，初步識(shí)別為干旱事件；②對(duì)于干旱歷時(shí)僅為1 個(gè)月的干旱事件，若指數(shù)值大于R2，則該干旱事件應(yīng)被剔除；③對(duì)于間隔時(shí)間為1 個(gè)月的2 個(gè)干旱事件，若間隔月份的干旱指數(shù)值小于R0，則將這2 個(gè)干旱事件合并為1 個(gè)干旱事件。 圖2 中干旱事件a、b分別識(shí)別為2 次獨(dú)立的干旱事件，但其間隔月份為1 個(gè)月，且該月干旱指數(shù)小于R0，因此將a、b及其間隔融合為1 次干旱事件進(jìn)行考慮，合并后，干旱歷時(shí)(烈度)為2 個(gè)干旱事件歷時(shí)(烈度)之和。 本文R0=0，R1=0.3，R2=0.5[22-24]，示意見(jiàn)圖2。

圖2 游程理論示意

2.3 干旱傳遞模型的構(gòu)建

Copula函數(shù)是由Sklar[25]在1959 年提出，已經(jīng)廣泛地應(yīng)用在氣象水文領(lǐng)域。 該理論的原理相對(duì)簡(jiǎn)單，認(rèn)為任意一個(gè)多維度聯(lián)合分布函數(shù)可以由一個(gè)Copula函數(shù)和多個(gè)邊緣分布函數(shù)構(gòu)成，并且允許多個(gè)邊緣分布函數(shù)的分布概率存在差異，可以用于描述任意變量間的聯(lián)合分布情況。 因此，擬采用二變量Copula函數(shù)構(gòu)建氣象-水文干旱傳遞模型。

二變量Copula模擬主要包括邊緣分布函數(shù)擬合及優(yōu)選、Copula函數(shù)擬合及優(yōu)選、聯(lián)合概率或條件概率分析等4 個(gè)過(guò)程。 以歷時(shí)為例，氣象-水文干旱傳遞模型的構(gòu)建方法為:①選擇伽瑪分布、對(duì)數(shù)正態(tài)分布、廣義正態(tài)分布、皮爾遜三型等單變量分布模型，對(duì)氣象干旱歷時(shí)和水文干旱歷時(shí)進(jìn)行邊緣分布擬合，確定各自最優(yōu)邊緣分布函數(shù)；②分析氣象干旱歷時(shí)和水文干旱歷時(shí)的相關(guān)性，若其相關(guān)性較強(qiáng)，符合構(gòu)建Copula函數(shù)的條件； ③選擇Gaussian、Gumbel、Frank 和Clayton 等常用二變量Copula函數(shù)，模擬氣象干旱歷時(shí)和水文干旱歷時(shí)聯(lián)合分布，選擇最優(yōu)Copula函數(shù)；④計(jì)算中旱、重旱、極旱3 個(gè)等級(jí)的水文干旱發(fā)生的氣象干旱條件概率，基于條件概率確定該等級(jí)干旱歷時(shí)傳遞閾值。

2.3.1 擬合邊緣分布

擬合邊緣分布是構(gòu)建Copula函數(shù)的首要任務(wù)，目前常見(jiàn)的4 種單變量函數(shù)是伽瑪分布、對(duì)數(shù)正態(tài)分布、廣義正態(tài)分布、皮爾遜三型分布。 本文設(shè)氣象干旱歷時(shí)和烈度的邊緣分布函數(shù)分別為FM(d)、FM(s)，水文干旱歷時(shí)和烈度的邊緣分布函數(shù)分別為FH(d)、FH(s)，采用上述邊緣分布函數(shù)對(duì)氣象干旱、水文干旱的歷時(shí)、烈度進(jìn)行擬合，用Anderson-Darling test(A-D)檢驗(yàn)法[26]進(jìn)行擬合優(yōu)度的檢驗(yàn)，認(rèn)為處于95%置信區(qū)間的統(tǒng)計(jì)量符合檢驗(yàn)要求，即通過(guò)檢驗(yàn)。 同時(shí)，采用Kendall秩相關(guān)系數(shù)、Pearson相關(guān)系數(shù)和Spearman 秩相關(guān)系數(shù)等方法表征干旱歷時(shí)和干旱烈度之間的相關(guān)性。

2.3.2 Copula擬合

Copula函數(shù)是連接單變量邊緣分布，形成在[0，1]區(qū)間上服從均勻分布的多元函數(shù)。 由Copula函數(shù)的定義可知，氣象干旱歷時(shí)FM(d)和水文干旱歷時(shí)FH(d)的聯(lián)合分布函數(shù)為:

同理可以計(jì)算氣象干旱烈度FM(s)和水文干旱烈度FH(s)的聯(lián)合Copula函數(shù)。 文中選取常用的理論二變量Copula函數(shù)Gaussian-copula和Claytoncopula、Frank-copula、Gumbel-copula[27]。 計(jì)算理論Copula值與經(jīng)驗(yàn)Copula值之間的均方根誤差RMSE[28]評(píng)價(jià)Copula函數(shù)的擬合優(yōu)度。

式中 N——樣本容量；Pc(i)——實(shí)測(cè)概率值；P0(i)——Copula函數(shù)計(jì)算得到的概率值。

2.3.3 條件概率

假設(shè)X、Y為2 個(gè)隨機(jī)變量，則X≥u 條件下Y≥v的條件概率可表示為:

式中 x、y——隨機(jī)變量X、Y的邊緣分布；x(u)、y(v)——X≤u 和Y≤v時(shí)的累積概率；C(x(u)，y(v))——X≤u 和Y≤v時(shí)的聯(lián)合概率。

本文X為氣象干旱特征，Y為不同等級(jí)的水文干旱特征。 干旱通常分為輕旱、中旱、重旱和極旱4個(gè)等級(jí)[29]，表1 所示，輕旱、中旱、重旱和極旱等級(jí)干旱事件的累積概率區(qū)間為(0，0.5)、[0.5，0.75)、[0.75，0.9)和[0.9，1)[30-31]。 構(gòu)造了氣象干旱特征和水文干旱特征的聯(lián)合分布模型后，根據(jù)式(6)計(jì)算氣象干旱特征大于某一個(gè)值時(shí)水文干旱等級(jí)超過(guò)中旱、重旱、極旱的條件概率。 隨著氣象干旱條件的增加，條件概率逐漸接近1。 本文選取干旱傳播閾值的置信水平為0.95[22]，即當(dāng)條件概率等于0.95 時(shí)，對(duì)應(yīng)的氣象干旱特征值被視為該等級(jí)水文干旱的干旱特征傳遞閾值，干旱等級(jí)劃分見(jiàn)表1。

表1 干旱等級(jí)劃分

3 結(jié)果分析

3.1 干旱指數(shù)特征

東江流域月尺度的SPI和SRI見(jiàn)圖3，大部分年份SPI值比SRI小。 SPI的最小值為-1.45，出現(xiàn)在1979 年10 月，最大值為2.55，出現(xiàn)在1966 年6 月。SRI的最小值為-1.60，出現(xiàn)在1963 年5 月，最大值為2.52，出現(xiàn)在1959 年6 月。 SPI的最大值和最小值均比SRI大。

圖3 東江流域月尺度的SPI和SRI

氣象干旱和水文干旱均表現(xiàn)為春旱(1—3 月)和冬旱(10—12 月)。 但氣象干旱一般出現(xiàn)于10 月份，于次年2 月份結(jié)束，平均持續(xù)時(shí)間為4 個(gè)月。 歷時(shí)最長(zhǎng)的氣象干旱出現(xiàn)在1962 年11 月至1963 年5月、1963 年10 月至1964 年4 月、2001 年10 月至2002 年4 月，歷時(shí)均為7 個(gè)月。 水文干旱一般出現(xiàn)于10 月份，于翌年3 月份結(jié)束，平均持續(xù)時(shí)間為5個(gè)月。 歷時(shí)最長(zhǎng)的水文干旱出現(xiàn)在1962 年11 月至1964 年4 月，歷時(shí)18 個(gè)月。 同時(shí)，對(duì)于氣象干旱，輕旱等級(jí)所占比例由55.5%下降至37.2%，中旱等級(jí)所占比例由44.6%上升至62.8%，進(jìn)入21 世紀(jì)后輕旱所占比例上升至47.5%，中旱等級(jí)所占比例減少至52.5%，與氣象干旱相比水文干旱變化規(guī)律不明顯。

3.2 干旱事件及其特征

計(jì)算干旱指數(shù)后，根據(jù)游程理論提取干旱事件。東江流域發(fā)生氣象干旱事件62 次，發(fā)生水文干旱事件54 次。 對(duì)時(shí)間上存在部分或全部重疊的干旱事件進(jìn)行匹配。 例如1956 年9 月至1957 年1 月發(fā)生歷時(shí)為5 個(gè)月的氣象干旱，1956 年9 月至1957 年2月發(fā)生歷時(shí)為6 個(gè)月的水文干旱，二者在發(fā)生時(shí)間上部分重疊，可匹配。 匹配后有干旱事件45 次。 匹配后干旱特征箱型見(jiàn)圖4，統(tǒng)計(jì)參數(shù)見(jiàn)表2。 氣象干旱歷時(shí)均值為4.8 個(gè)月，最小值為1 個(gè)月，最大值為7 個(gè)月；烈度均值為3.24，最小值為0.91，最大值為4.93；水文干旱歷時(shí)均值為5.42 個(gè)月，最小值為1個(gè)月，最大值為8 個(gè)月；烈度均值為3.11，最小值為0.35，最大值為5.42。 水文干旱相對(duì)氣象干旱平均歷時(shí)長(zhǎng)，烈度小。

圖4 干旱特征箱型

表2 干旱特征的統(tǒng)計(jì)參數(shù)

3.3 氣象-水文干旱聯(lián)合分布擬合

識(shí)別干旱事件后，采用伽瑪分布、對(duì)數(shù)正態(tài)分布、廣義正態(tài)分布、皮爾遜三型分布對(duì)匹配后的氣象干旱和水文干旱的歷時(shí)和烈度分別進(jìn)行擬合，在通過(guò)顯著性檢驗(yàn)的前提下，結(jié)合歷時(shí)和烈度的累積概率曲線圖，根據(jù)AIC最小原則，選擇最優(yōu)擬合函數(shù)。各個(gè)邊緣分布函數(shù)對(duì)歷時(shí)和烈度的擬合效果見(jiàn)圖5，擬合參數(shù)見(jiàn)表3。 氣象干旱歷時(shí)和烈度的最優(yōu)擬合函數(shù)均為伽瑪分布，水文干旱歷時(shí)的最優(yōu)邊緣分布函數(shù)為伽瑪分布，水文干旱烈度的邊緣分布函數(shù)為廣義正態(tài)分布。 采用Kendall秩相關(guān)系數(shù)、Pearson 相關(guān)系數(shù)和Spearman 秩相關(guān)系數(shù)計(jì)算氣象干旱歷時(shí)(烈度)和水文干旱歷時(shí)(烈度)之間的相關(guān)性，結(jié)果見(jiàn)表4，變量呈明顯的正相關(guān)關(guān)系，其中，氣象干旱烈度與水文干旱烈度之間的相關(guān)性較高。

圖5 邊緣分布函數(shù)擬合效果

表3 邊緣分布函數(shù)統(tǒng)計(jì)檢驗(yàn)P值

表4 氣象干旱-水文干旱特征相關(guān)系數(shù)

建立氣象干旱歷時(shí)(或烈度)和水文干旱歷時(shí)(或烈度)的二變量Copula函數(shù)，并計(jì)算理論Copula值和經(jīng)驗(yàn)Copula值的RMSE驗(yàn)證Copula函數(shù)的擬合效果。 表5 為評(píng)價(jià)指標(biāo)RMSE的值，其中RMSE越小，說(shuō)明Copula函數(shù)擬合效果越好。 對(duì)于烈度，Gaussian-copula的RMSE值最小，為0.047 2；對(duì)于歷時(shí)，Gumbel-copula的RMSE值最小，為0.094 6。 說(shuō)明Copula函數(shù)擬合效果滿足預(yù)期要求，且滿足程度較高。 由于不同Copula函數(shù)的結(jié)構(gòu)不同，會(huì)導(dǎo)致Copula函數(shù)的選擇影響分析結(jié)果。 因此，本文統(tǒng)一采用Gumbel-copula對(duì)干旱特征進(jìn)行擬合。

表5 Copula 函數(shù)擬合效果評(píng)價(jià)指標(biāo)

3.4 氣象干旱向水文干旱傳遞分析

本文構(gòu)建Copula函數(shù)后，計(jì)算了給定氣象干旱歷時(shí)(烈度)超過(guò)某一值條件下水文干旱等級(jí)大于中旱(累積概率大于0.5)、重旱(累積概率大于0.75)和極旱(累積概率大于0.9)的概率分布。 不同氣象干旱特征下，中旱、重旱和極旱等級(jí)水文干旱發(fā)生的條件概率見(jiàn)圖6。 結(jié)果表明，相同氣象干旱特征下，水文干旱的條件概率隨著干旱等級(jí)的增加而減小。 條件概率曲線呈“S”形，在條件概率接近0.95 時(shí)，隨著氣象干旱特征的增加，條件概率曲線趨于平緩。 截取條件概率為0.95 時(shí)的干旱特征值作為干旱傳遞的閾值。 東江流域中旱、重旱和極旱等級(jí)的水文干旱歷時(shí)和烈度的閾值為6.9、8.3、9.4個(gè)月和4.6、5.7、6.6。 隨著干旱等級(jí)的增加，干旱的閾值也相應(yīng)增加。

圖6 氣象干旱向水文干旱傳遞的條件概率

4 結(jié)論

計(jì)算東江流域的SPI指數(shù)和SRI指數(shù)，根據(jù)游程理論提取干旱歷時(shí)和烈度，對(duì)干旱特征進(jìn)行邊緣分布擬合，建立氣象干旱歷時(shí)(烈度)和水文干旱歷時(shí)(烈度)的二元Copula模型，計(jì)算不同氣象干旱條件下中旱、中旱、極旱水文干旱發(fā)生的條件概率，截取條件概率為0.95 時(shí)的氣象干旱特征值作為干旱傳遞閾值。 主要結(jié)論如下。

a)SPI和SRI分別體現(xiàn)了氣象干旱和水文干旱嚴(yán)重程度，根據(jù)游程理論識(shí)別了氣象干旱62 次，水文干旱54 次，歷時(shí)最長(zhǎng)的氣象干旱為7 個(gè)月，出現(xiàn)在1962 年11 月至1963 年5 月、1963 年10 月至1964 年4 月、2001 年10 月至2002 年4 月，歷時(shí)最長(zhǎng)的水文干旱為18 個(gè)月，出現(xiàn)在1962 年11 月至1964年4 月，匹配后干旱事件45 次。

b)通過(guò)A-D檢驗(yàn)方法對(duì)邊緣分布擬合結(jié)果進(jìn)行檢驗(yàn)，氣象干旱歷時(shí)和烈度的最優(yōu)擬合函數(shù)均為伽馬分布，水文干旱歷時(shí)和烈度的最優(yōu)擬合函數(shù)分別為伽瑪分布、廣義正態(tài)分布。 通過(guò)RMSE結(jié)果對(duì)比4 種Copula函數(shù)的擬合效果，Gumbel-copula函數(shù)的擬合效果最優(yōu)，則選擇Gumbel-copula函數(shù)構(gòu)建干旱傳遞模型。

c)氣象-水文干旱特征聯(lián)合分布條件概率取0.95 時(shí)，東江流域中旱、重旱和極旱的干旱歷時(shí)傳遞閾值分別為6.9、8.3、9.4 個(gè)月，干旱烈度傳遞閾值分別為4.6、5.7、6.6。