東江流域氣象干旱向水文干旱傳遞閾值研究

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

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

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

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

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

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

東江為珠江流域三大水系之一，流域控制站博羅站以上的集水面積為25 325 km2。 東江流域屬于亞熱帶季風氣候，流域多年平均降水量為1 500 ～2 400 mm，降雨時間上分配不均，4—9 月為雨季，降水量達全年降水量的80%，10 月到翌年3 月為旱季，降水量僅為全年降水量的20%。 東江擔負著為珠三角地區(qū)4 000 多萬人口提供生產、生活、生態(tài)用水的重任，因此研究東江流域干旱情況對社會和經濟發(fā)展有著重要意義。 研究數(shù)據(jù)為中國國家氣象局氣象信息中心提供的東江流域內34 個主要氣象站1956 年4 月至2009 年3 月的月降水數(shù)據(jù)和博羅水文站1956 年4 月至2009 年3 月的實測月徑流數(shù)據(jù)，并根據(jù)流域三大水庫(新豐江水庫、楓樹壩水庫和白盆珠水庫)的月徑流調蓄量，還原得到天然徑流數(shù)據(jù)。 站點分布見圖1。

圖1 東江流域及站點分布

2 研究方法

2.1 干旱指數(shù)計算

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

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

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

對累積概率進行標準化轉換，有:

因此，干旱指數(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ù)的計算理論為基礎，以徑流值代替原SPI計算過程的降水值，即可得到標準化徑流指數(shù)SRI[17]。

2.2 干旱事件識別

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

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

圖2 游程理論示意

2.3 干旱傳遞模型的構建

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

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

2.3.1 擬合邊緣分布

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

2.3.2 Copula擬合

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

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

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

2.3.3 條件概率

假設X、Y為2 個隨機變量，則X≥u 條件下Y≥v的條件概率可表示為:

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

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

表1 干旱等級劃分

3 結果分析

3.1 干旱指數(shù)特征

東江流域月尺度的SPI和SRI見圖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 月份結束，平均持續(xù)時間為4 個月。 歷時最長的氣象干旱出現(xiàn)在1962 年11 月至1963 年5月、1963 年10 月至1964 年4 月、2001 年10 月至2002 年4 月，歷時均為7 個月。 水文干旱一般出現(xiàn)于10 月份，于翌年3 月份結束，平均持續(xù)時間為5個月。 歷時最長的水文干旱出現(xiàn)在1962 年11 月至1964 年4 月，歷時18 個月。 同時，對于氣象干旱，輕旱等級所占比例由55.5%下降至37.2%，中旱等級所占比例由44.6%上升至62.8%，進入21 世紀后輕旱所占比例上升至47.5%，中旱等級所占比例減少至52.5%，與氣象干旱相比水文干旱變化規(guī)律不明顯。

3.2 干旱事件及其特征

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

圖4 干旱特征箱型

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

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

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

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

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

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

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

表5 Copula 函數(shù)擬合效果評價指標

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

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

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

4 結論

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

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

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

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