基于Copula函數和SPEI的廣東北江流域干旱頻率分析*

劉占明，陳子燊，黃 強，劉曾美

(1.中山大學水資源與環境系，廣東 廣州510275；2.華南理工大學水利水電系，廣東 廣州510640)

干旱是世界上最常見的自然災害，被認為是全球最嚴重的自然災害之一[1]。在我國，干旱被列為氣象災害之首[2]。雖然廣東位于華南沿海，降水較多，但在季風、地質條件等多種因素的綜合影響下，也飽受干旱困擾。根據統計[3]，廣東省約24％的耕地經常遭遇干旱，嚴重的年份耕地受旱面積超過100萬hm2(約占耕地總面積的50％)。

北江起源于南嶺山系，是珠江水系第二大支流，流域面積約4.7萬 km2，其中廣東境內約4.3萬km2。由于干旱事件具有多種屬性，如Mishra和Singh[4]采用歷時、烈度、嚴重程度等描述干旱，并且這些特征屬性間具有一定關聯，適宜進行多變量分析。從過去的研究來看，涉及該區域干旱的研究已有較多文獻[5-7]，但這些研究多是基于某一評估指標(主要是標準化降水指數SPI)研究干旱某一變量的時空變化特征，較少涉及多變量的概率分析。有研究[8]指出SPI沒有考慮溫度變化造成的水分支出，對干旱的評估結果還需進一步檢驗。北江流域屬中亞熱帶向南亞熱帶過渡區域，溫度變化造成的水分支出對于干旱的影響是不容忽視的。近年來，已有學者提出氣溫對于干旱的重要性并對蒸散發進行了研究與估算[9]。Vicente-Serrano等[10]提出了標準化降水蒸散指數SPEI(the standardized precipitation evapotranspiration index，簡稱SPEI)。SPEI綜合了帕默爾指數(PDSI)對蒸散的響應及SPI多時間尺度的優點，基于概率計算，具有空間可比性的優點，目前已成為研究干旱新的理想指標[11-12]。劉占明等[13]應用7種干旱指標對廣東北江流域進行評估，并結合歷史干旱資料對比分析，發現SPEI指標的評估結果最能反映該區域的實際情況。眾所周知，在氣象水文學領域，極端事件的各屬性之間普遍存在著相關性，而近10年來在氣象水文學領域引入的Copula函數可以描述這種相關性，并且為多變量聯合概率分布的研究提供了理論基礎和技術手段[14]。

本文選取廣東北江流域分布較均勻、數據序列較長且沒有發生站點變更的18個國家基本氣象站(圖1)近50a日均溫、日降水數據，數據序列大多起始于1957-1962年，均終止于2007年，數據資料來自國家氣象局氣象信息中心，均經過了嚴格的質量控制及三性檢查，具有良好的完整性。根據計算出的SPEI值，通過游程理論[15]確定干旱變量，基于概率函數對區域干旱單變量概率分布和兩變量聯合分布及其帶來的干旱風險進行研究，以期為分析區域水量平衡及防災減害的風險管理提供科學依據。

圖1 廣東北江流域范圍及站點位置

1 研究方法

1.1 標準化蒸散指數(SPEI)

基于水量平衡原理提出的SPEI，首先利用Thornthwaite[16]方法計算潛在月蒸散量，進而以潛在蒸散量與降水量差值的概率分布來表征一個區域凈降水量偏離正常的程度。計算方法如下：

1)計算地表月潛在蒸散量

EP=16K(10T/I)a

(1)

式中，Ep為潛在蒸散發量(mm)；T為月平均氣溫(℃)；I為年熱指數，等于一年中各月指數i=(T/5)1.514的累加值；a為經驗指數，由與I的函數關系導出：a=6.75×10-7I3-7.71×10-5I2+1.79×10-2I+0.49；修正系數K=(N/12)(M/30)，其中，M為每月的天數，N為最大日照時數，N=(24/π)ω，ω為日落時太陽時角，由緯度φ和太陽傾角δ決定：ω=arcos(-tanφtanδ)，δ=0.409 3sin(2πJ/365-1.405)，J為日序數。

2)計算逐月降水與蒸散的差值

Di=Pi-EP(i)

(2)

式中，Di為第i個月的凈降水量，mm。不同時間尺度的Di為：

(3)

式中k為月時間尺度。

3)計算Di的累積概率

采用3參數對數邏輯斯特分布(Log-Logistic)擬合序列Di，其分布函數為：

F(x)= 1/(1 + (α/(x-γ))β)

(4)

式中α、β和γ分別為尺度、形狀和位置參數，由概率權重矩法估計：

α=(w0-2w1)β/[Г(1+1/β)Г(1-1/β)]，

β=(2w1-w0)/(6w1-6w0w2)，γ=w0-Г(1+1/β)Г(1-1/β)，wk為k階樣本概率權重矩，k=0，1，2…。

4)累積概率標準正態化

凈降水量服從的Log-Logistic分布為偏態分布，將偏態分布累積概率值轉換為標準正態分布：

z=W-(c0+c1W+c2W2)/

(1+d1W+d2W2+d3W3)

(5)

式中，z為標準正態化后的SPEI值；c0=2.515517；c1=0.802853；c2=0.010328；d1=1.432788；d2=0.189269；d3=0.001308；當凈降水量累積概率F≤0.5時，W= (-2ln(1-F))1/2，當F>0.5時，W= (-2ln(F))1/2，相應的z值取相反數。

表1 SPEI干旱等級劃分

由于SPEI 通過概率密度函數求解累積概率，再將累積概率標準正態化而得，消除了降水和氣溫的時空分布差異，在不同區域和時段均能有效地反映旱澇狀況。與SPI一樣，SPEI是基于概率計算得到的，因此SPI的干旱等級劃分標準SPEI同樣適用，根據中華人民共和國《氣象干旱等級》[16]國家標準來劃分干旱等級，如表1所示。

SPEI可根據不同尺度進行計算，一般有1、3、6、12個月時間尺度。由于北江流域一年中降水主要集中在汛期(4-9月)，汛期又是溫度較高的夏半年，蒸發旺盛，非汛期(10月-次年3月)降水較少，蒸發較弱，時間尺度大致為6個月，因此為較好的識別干旱情況，本文選用6個月時間尺度的SPEI。

1.2 干旱識別與干旱變量的定義

圖2 游程理論識別干旱變量示意圖

1.3 干旱變量邊緣分布

一般來說，對于干旱歷時而言，如果作為離散型分布，可用幾何分布擬合[18]，如果作為連續型分布，常用指數分布進行擬合[17]，結合當前國內外的研究情況，本文采用指數分布對干旱歷時進行擬合。在干旱烈度和烈度峰值的擬合上，本文選用氣象水文學領域常用的指數分布(兩參數，即2P，下同)、標準耿貝爾分布(2P)、伽瑪分布(皮爾遜 Ⅲ型，3P)、廣義極值分布(3P)、廣義邏輯斯特分布(3P)、對數邏輯斯特分布(3P)、對數正態分布(3P)、廣義帕累托分布(3P)、韋克比分布(5P)、維布爾分布(3P)、銳利分布(2P)共11種分布函數對18站進行擬合。統一采用線性矩法估計參數[19]，通過Kolmogorov-Smirnov(K-S)法進行優度檢驗[20]，并確定出最佳的分布函數，進而運用Copula函數構建聯合分布。

1.4 Copula函數定義

1.5 Copula函數的選擇與參數估計

本文根據Genest和Rivest[21]提出的理論估計值-經驗估計值關系圖法來評價和選擇Copula。采用Kendall秩相關系數τ[14]度量連接函數的相關性，并利用相關系數τ計算出二維Archimedean Copula函數的參數值。

Shiau和Shen[22]推導出了干旱變量聯合重現期上(T0)和同現重現期(Ta)的計算公式：

T0(x,y) =E(L)/P[X>xorY>y] =

E(L)/[1-C(FX(x),FY(y))]

(6)

Ta(x,y) =E(L)/P[X>x,Y>y]

=E(L)/[1-FX(x)-FY(y)+C(FX(x),FY(y))]

(7)

式中，E(L)為干旱間隔期望值(即非干旱歷時與干旱歷時期望值之和)。變量X和Y的重現期(即邊緣重現期)計算公式為：

T(x) =E(L)/[1-FX(x)]；

T(y) =E(L)/[1-FY(y)]

(8)

以上計算過程均采用Matlab軟件編程實現，部分代碼由http://code.google.com/p/copula/網站免費獲得。

2 結果與分析

2.1 邊緣分布函數結果分析

2.1.1 概率分布函數選擇 研究結果表明(表2所示)，干旱烈度的擬合中兩參數的銳利分布、指數分布、標準耿貝爾分布擬合較差，分別有83.3％、55.6％、16.7％的站沒有通過檢驗(0.05置信水平，下同)，其它分布除對數正態分布(3P)有1個站(仁化)沒有通過檢驗外，其余都通過了檢驗，根據K-S D統計值越小擬合越優原則，可以看出維布爾分布(3P)、韋克比分布(5P)、伽瑪分布(3P)、對數正態分布(3P)和廣義帕累托分布(3P)在大部分站點擬合較好，其中維布爾分布(3P)在大多數站點擬合最優。烈度峰值的擬合中(由于篇幅所限，不再以表格形式列出)，18站在所有分布函數中都能通過0.05臨界值檢驗，其中維布爾分布(3P)和韋克比分布(5P)在大部分站點擬合較好，而廣義帕累托分布(3P)在大多數站點擬合最優。

表2 干旱烈度11種概率分布的K-S D統計量1)

1)加粗的為擬合較好的前3種分布函數，標*的為擬合最好的概率分布(也是本文采用的概率分布)。

2.1.2 干旱變量10a重現期設計值空間分布 表3為根據上述擬合最好的概率分布函數推斷出的干旱烈度(S)、烈度峰值(P)及干旱歷時(D，指數分布)5、10、20、50、100 a一遇設計值(只列出了韶關站數據)；這里主要分析干旱歷時、烈度和烈度峰值10年一遇設計值的空間分布情況(圖3所示，采用樣條函數法[23]插值，該方法的優點是能夠準確的通過樣點擬合出連續光滑的表面)。

表3 韶關站各重現期對應的干旱變量設計值及聯合分布的重現期

圖3 干旱歷時(a)、烈度(b)及烈度峰值(c)10年一遇設計值空間分布

流域干旱歷時10 a重現期設計值(圖3a)最大值為始興9.3月，最小值為英德7.7月，其余各站均介于8-9月之間，流域18站平均為8.5月；總體來看流域內部區域差異并不明顯，但平均8.5月的干旱歷時說明流域面臨干旱歷時較長的風險較大。

干旱烈度10 a重現期設計值(圖3b)最大值(10.13)位于陽山，其余各站均介于8.5～10.0之間，平均值為9.19；總體來看，流域中部偏西及北部邊緣(陽山、樂昌、仁化等地區)為高值區域，東部的英德、翁源地區為低值區域。

烈度峰值10 a重現期設計值(圖3c)最小值為1.62(陽山站)，連山、英德、始興、懷集、廣寧、三水、清遠7站(占流域18站38.89％)設計值在2.0以上，18站平均值為1.96；從中可以看出流域18站10年一遇設計值均達到重旱以上級別，平均值接近特旱級別，38.89％的站達到特旱(等級劃分標準見文獻[13，16])級別；空間分布趨勢表現為西南部邊緣及東部為高值區域，中部偏西及北部邊緣為低值區域。

2.2 二維聯合Copula函數結果分析

2.2.1 Copula函數的選擇與參數估計結果分析 計算發現18站干旱歷時、干旱烈度、烈度峰值兩兩之間的Kendall秩相關系數τ分別介于0.789～0.903、0.519～0.706、0.673～0.848之間，說明干旱變量之間具有較高的相關性。由于AMH Copula函數不適用于高的正或負相關性，只有Kendall系數-0.1817<τ<1/3時才適用，而本研究中的隨機變量Kendall系數均超過1/3，故不適用于AMH Copula函數。

由Genest-Rivest檢驗法[21]計算的Kc(t)-Ke(t)關系圖(圖4只列出了韶關站Kc(t)-Ke(t)關系圖)發現，18站干旱歷時-烈度、歷時-烈度峰值均表現為GH Copula函數圖上的點最接近45°對角線，即GH Copula擬合最優；18站烈度-烈度峰值Kc(t)-Ke(t)關系圖中，GH Copula、Clayton Copula函數圖上的點均分布在45°對角線附近，其中8站表現為GH Copula函數圖上的點最接近45°對角線，10站表現為Clayton Copula函數圖上的點最接近45°對角線。在干旱風險分析中，人們更關心的是隨機變量的尾部相關性，也就是當一個隨機變量取較大的值或者較小的值時，它對另一個隨機變量的取值是否有影響。根據干旱的理論知識可知，干旱歷時的長短是干旱烈度及烈度峰值的重要原因，由于干旱具有連續性的累積效應，一般來說干旱歷時越長，烈度越大，旱情越嚴重，對應的烈度峰值也會越大；結合干旱的風險評估，人們更關注干旱歷時、烈度、烈度峰值之間的上尾相關性，因此在烈度-烈度峰值GH Copula與Clayton Copula的對比中本文選取具有上尾相關性的GH Copula作為聯接函數。然后采用相關性指標法估計聯接函數的參數θ，即運用Kendall秩相關系數τ與GH Copula參數θ之間的關系(τ=1-1/θ)計算。

圖4 韶關站Genest-Rivest方法檢驗結果

2.2.2 Copula兩變量重現期和風險分析 表3列出了韶關站5、10、20、50、100年重現期對應的干旱變量二維聯合分布(D-S、D-P、S-P)的聯合重現期(T0)和同現重現期(Ta)，本文主要分析流域50和100 a的情況(圖5，采用樣條函數插值法[23])。

從圖5可以看出，流域D-S50、100 aT0(圖5a1、a2)及Ta(圖5 b1、b2)、D-P50、100 aT0(圖5c1、c2)、S-P50、100 aT0(圖5e1、e2)具有比較相似的空間分布趨勢，均主要表現為西部(連南、連山、懷集、廣寧一帶)、南部(三水、清遠地區)及北部的始興地區為相對低值區域，中北部的陽山、樂昌、仁化及東部的翁源地區為相對高值區域。D-P50、100 aTa(圖5d1、d2)、S-P50、100 aTa(圖5f1、f2)具有比較相似的空間分布趨勢，雖然低值區域仍然為西部、南部及北部的始興地區，但其高值區域范圍明顯擴大，北部(樂昌、南雄地區)、東部(翁源、佛岡地區)及中部的乳源、韶關、英德地區均為高值區域。可以看出，流域西部(連南、連山、懷集、廣寧一帶)、南部(三水、清遠地區)及北部的始興地區出現其中之一(即干旱歷時長或烈度大或烈度峰值大)或其中兩者同時出現的概率均較大，從而具有較大的干旱風險；相對而言，流域中北部(陽山、韶關、樂昌等地)及東部的翁源等地干旱風險較小。

從流域18站平均值來看，D-S50、100 aT0平均值分別為42.3、84.5 a，Ta平均值分別為52.1和104.0 a；D-P50、100 aT0平均值分別為35.3和70.4 a，Ta平均值分別為69.2和138.9 a；S-P50、100年T0平均值分別為39.0、78.4 a，Ta平均值分別為57.7、115.9 a。從最大值最小值來看，D-S、D-P及S-P50、100 aT0最小值均位于懷集，最大值均位于陽山；D-S、D-P及S-P50、100 aTa最小值均位于懷集，最大值均位于翁源。從中可以看出，整個流域出現干旱歷時長或烈度峰值大的概率較大，出現干旱烈度大或烈度峰值大的概率也較大，其中懷集出現以上各種干旱情況的概率最大。

圖5 廣東北江流域干旱變量聯合分布重現期空間分布

3 結 論

本文基于SPEI指標對廣東北江流域18站近50 a干旱單變量概率分布和兩變量Copula聯合分布進行研究，主要得到以下結論：

1) 干旱烈度的擬合中兩參數的銳利分布、指數分布擬合較差，韋克比分布(5P)、伽瑪分布(皮爾遜Ⅲ，3P)、對數正態分布(3P)和廣義帕累托分布(3P)在大部分站點擬合較好，維布爾分布(3P)在大部分站點擬合最優；烈度峰值的擬合中，維布爾分布(3P)和韋克比分布(5P)在大部分站點擬合較好，廣義帕累托分布(3P)在大部分站點擬合最優。

2) 干旱歷時十年一遇設計值區域差異不明顯，平均值為8.5月，說明整個流域面臨干旱歷時較長的風險較大。流域18站烈度峰值十年一遇設計值均達到重旱以上級別，平均值接近特旱級別，近40％的站達到特旱級別；由于烈度峰值是干旱程度的極端表現，以上結論說明流域在極端干旱上形勢嚴峻。

3) 從流域50、100年重現期對應的干旱變量二維聯合分布的聯合重現期和同現重現期來看，流域西部(連南、連山、懷集、廣寧一帶)、南部(三水、清遠地區)及北部的始興地區出現其中之一(即干旱歷時長或烈度大或烈度峰值大)或其中兩者同時出現的概率均較大，從而具有較大的干旱風險；相對而言，流域中北部(陽山、韶關、樂昌等地)及東部的翁源等地干旱風險較小。

參考文獻：

[1]WIHITE DA.Drought as a Natural Hazard:Concepts and Definitions[M]//WIHITE D A,Ed.Drought:A Global Assessment.Routledge,2000:3-18.

[2]國家科學技術委員會．中國科學技術藍皮書第5號《氣候》[M]．北京：科學文獻出版社，1990.

[3]中國氣象災害大典·廣東卷[M].北京:氣象出版社,2006:6,204-206.

[4]MISHRA ASHOK K,SINGH VIJAY P.A review of drought concepts[J].Journal of Hydrology,2010,391:202-216.

[5]黃晚華，楊曉光，李茂松,等．基于標準化降水指數的中國南方季節性干旱近58a演變特征[J]．農業工程學報，2010，26(7):50-59.

[6]陳子燊，黃 強，劉曾美．1962-2007年廣東干濕時空變化特征分析[J]．水科學進展，2013，24(4):469-476.

[7]江 濤，楊 奇，張 強,等．廣東省干旱災害空間分布特征[J]．湖泊科學，2012，24(1):156-160.

[8]McKEE TB，DOESKEN NJ，KLEIST J．The relationship of drought frequency and duration to time scales[C]//The 8th Conference on Applied Climatology，Anaheim，CA，1993:179-184.

[9]HU Q,WILLSON G D.Effect of temperature anomalies on the Palmer drought severity index in the central United States[J].International Journal of Climatology,2000,20,1899-1911.

[10]VICENTE-SERRANO SM,BEGUERíA S,LóPEZ-MORENO JI.A multiscalar drought index sensitive to global warming:The standardized precipitation evapotranspiration index[J].Journal of Climate,2010,23(7):1696-1718.

[11]CAI W，COWAN T．Evidence of impacts from rising temperature on inflows to the Murray-darling basin[J]．Geophysical Research Letters，2008，35:L07701．[doi:10．1029/2008GL033390].

[12]LORENZO-LACRUZ J，VICENTE-SERRANO SM，LPEZ-MORENO J I，et al．The impact of droughts and water management on various hydrological systems in the headwaters of Tagus River (central Spain) [J]．Journal of Hydrology，2010，386:13-26.

[13]劉占明，陳子燊，黃 強，等．7 種干旱評估指標在廣東北江流域應用中的對比分析[J]．資源科學，2013，35(5):1007-1015.

[14]NELSEN RB.An Introduction to Copulas [M].New York:Springer,2006.

[15]EMIR Z,ATILA S.A method of streamflow drought analysis [J].Water Resources Research,1987,23(1) :156-168.

[16]氣象干旱等級 [S].GB/T20481-2006.

[17]SHIAU JT.Fitting drought duration and severity with two-dimensional Copulas[J].Water Resources Management,2006,20:795-815.

[18]MATHIER L,PERREAULT L,BOBE B.The use of geometric and gamma-related distributions for frequency analysis of water deficit[J].Stochastic Hydrology and Hydraulics,1992,6(4):239-254.

[19]HOSKING JRM.L-Moments:Analysis and estimation of distributions using linear combinations of order statistics[J].Journal of the Royal Statistical Society,1990,52(1):105-124.

[20]FRANK J,MASSE J.The Kolmogorov-Smirnov test for goodness of fit[J].Journal of the American Statistical Association,1951,46(253):68-78.

[21]GENEST C，RIVEST L．Statistical inference procedures for bivariate Archimedean copulas[J]．Journal of American Statistical Association，1993，88:1034-1043.

[22]SHIAU JT,SHEN HW.Recurrence analysis of hydrologic droughts of differing severity[J].Journal of Water Resources Planning and Management,2001,127(1):30-40.

[23]FRANRE R.Smooth interpolation of scattered path by local thin plate [J].Comp Maths with Appls Great Britain,1982,8(4):237-281.