顧佳帥，李瓊芳,2，牛銘媛，陳啟慧，和鵬飛，周正模，曾天山，杜 堯，宋 耘，韓幸燁

(1.河海大學水文水資源學院，南京 210098；2.江蘇省“世界水谷”與水生態(tài)文明協(xié)同創(chuàng)新中心，南京 210000)

0 引 言

特殊的地理位置以及氣候條件導致我國成為一個旱災頻發(fā)的國家,其中淮河上游流域地處我國南北氣候過渡帶，氣候條件復雜，旱災頻發(fā)[1]。因此，科學表征干旱事件對淮河上游流域的防旱抗旱具有重要的意義。干旱事件通常包含多方面的特征屬性，如干旱歷時、干旱烈度和干旱烈度峰值等[2]，但以往的干旱特征分析通常僅針對某單一干旱特征變量，然而單變量分析不能很好地表征干旱事件特征。相對單變量而言，多變量聯(lián)合分析充分考慮了干旱特征變量間的相互關系，因而分析結果更加全面[3]。傳統(tǒng)的多變量聯(lián)合分析方法對變量的邊緣分布有嚴格的要求，如只有具有相同邊緣分布多維變量才可采用多變量正態(tài)分布模型等經驗分布模型來擬合聯(lián)合分布[4]。基于Copula函數(shù)的多變量聯(lián)合分布模型不需要考慮變量邊緣分布的類型[2]，因此該模型在干旱等研究領域得到了廣泛應用。陸桂華等[5]以重慶市為研究區(qū)，利用Gumbel Copula函數(shù)構建干旱歷時與干旱烈度的二維聯(lián)合分布并計算重現(xiàn)期；馬明衛(wèi)、宋松柏[6]以西安站為例，利用橢圓形Copula函數(shù)構建干旱歷時、干旱烈度及烈度峰值三者的聯(lián)合分布并計算條件重現(xiàn)期；李計等[7]以烏魯木齊和石河子為研究區(qū)，利用Copula函數(shù)構建二維及三維干旱特征變量聯(lián)合分布模型并計算重現(xiàn)期；大量相關研究均表明Copula函數(shù)在構建多維干旱特征變量聯(lián)合分布函數(shù)方面具有優(yōu)勢。然而，國內外學者在淮河上游流域干旱特征分析的研究中極少運用Copula函數(shù)來構建干旱特征變量的聯(lián)合分布，因此，論文選取淮河上游息縣以上流域為研究區(qū)，計算月尺度長系列標準化降水指數(shù)(SPI)，利用三參數(shù)分布函數(shù)擬合干旱歷時、烈度以及烈度峰值三個干旱特征變量的邊緣分布，進而基于Copula函數(shù)分別構建其二維和三維聯(lián)合分布模型，推求多維變量的組合重現(xiàn)期并以此分析研究區(qū)的干旱特征，以期為流域防旱抗旱提供理論依據。

1 基于Copula函數(shù)的多維干旱變量聯(lián)合分布模型的構建

1.1 干旱特征變量的識別與確定

由于雨量站在流域上呈現(xiàn)非規(guī)則的離散分布，難以完全實際地反映降雨在流域的分布情況。因此研究采用泰森多邊形法推求流域的面雨量值,再由面雨量值進一步推求標準化降水指數(shù)。

標準化降水指數(shù)SPI計算方法較為簡單，能夠較好地反映干旱強度和干旱持續(xù)時間，在我國有很好的適用性[8]。干旱特征變量的定義和提取以游程理論[9]為基礎。干旱歷時D代表單次干旱過程的持續(xù)時間；干旱烈度S代表單次干旱過程持續(xù)時間內SPI的累加值；干旱烈度峰值P代表單次干旱過程中SPI的最小值。

識別的具體步驟為：確定SPI閾值R0、R1和R2，分別為-0.5、-1和0。若單月SPI值小于R1，則判定為一次干旱過程[如圖1(a)]，D為1個月，S為a的陰影面積，P在數(shù)值上等于S；若連續(xù)n月的SPI值均小于R0，則算做一次干旱過程[如圖1(e)]，D為n，S為n月SPI的和，P為SPI的最小值；若兩次干旱過程的時間間隔為一個月且該月的SPI值小于R2，則兩次干旱過程合并為一次干旱過程[如圖1(b)，(c)，(d)]，D為Db+1+Dd，S為Sb+Sd，P為SPI的最小值。

圖1 干旱特征變量的定義和提取Fig.1 Definition and extraction of drought characteristic variables

1.2 干旱特征變量邊緣分布函數(shù)的確定

采用對數(shù)正態(tài)分布、三參數(shù)伽馬分布以及三參數(shù)韋伯分布擬合干旱歷時、干旱烈度以及干旱烈度峰值的邊緣分布。經擬合優(yōu)度檢驗后，干旱歷時D采用三參數(shù)伽馬分布擬合，干旱烈度S以及干旱烈度峰值P的邊緣分布采用三參數(shù)韋伯分布擬合。

三參數(shù)伽馬分布的分布函數(shù)[10]為：

(1)

式中：a是形狀參數(shù)；b是尺度參數(shù)；c是位置參數(shù)。

三參數(shù)韋伯分布的分布函數(shù)[11]為：

(2)

式中：α是位置參數(shù)；β是形狀參數(shù)；η是尺度參數(shù)。

采用極大似然法估計參數(shù)，擬合優(yōu)度評價則采用Kolmogorov-Smimov(K-S)檢驗法[12]，當K-S檢驗統(tǒng)計量小于臨界值時，表示檢驗通過。

1.3 多維干旱變量聯(lián)合分布模型的建立

Copula函數(shù)是連接單變量邊緣分布，形成在[0，1]區(qū)間上服從均勻分布的多元函數(shù)[13]。設F1，F(xiàn)2，…，F(xiàn)n為連續(xù)隨機變量，則有唯一Copula函數(shù)C對任意x∈Rn滿足[14]：

F(x1,x2,…,xn)=P{X1≤x1,X2≤x2,…,Xn≤xn}=

C[F1(x1),F2(x2),…,Fn(xn)]

(3)

式中：x1,x2,…,xn為觀測樣本；F(x)為邊緣分布函數(shù)。

阿基米德型Copula函數(shù)結構相對簡單，計算方便，可以構造出多種形式、適應性強的多變量聯(lián)合分布函數(shù)，在水文領域的實際應用中較為廣泛。[15]本文采用阿基米德Copula函數(shù)中的四種：Frank、Clayton、Gumbel以及Joe Copula函數(shù)構建多維干旱變量的聯(lián)合分布。參數(shù)θ的估計采用半參數(shù)極大似然估計法[16]。擬合優(yōu)度的檢驗選用均方根誤差準則(RMSE)和赤池信息量準則(AIC)[17]， 以RMSE和AIC值最小作為擬合最優(yōu)的檢驗標準，公式為：

(4)

式中：MSE是均方誤差；n是樣本容量；i是樣本編號；Pei是經驗概率值；Pi是Copula函數(shù)聯(lián)合分布理論概率值；l是Copula函數(shù)的參數(shù)個數(shù)。

2 多維干旱變量聯(lián)合分布模型在淮河上游的應用

2.1 研究區(qū)概況

選取淮河上游息縣以上流域為研究區(qū)，流域范圍：東經113°15′～114°46′，北緯31°31′～ 32°43′，流域面積10 190 km2。該流域地處我國半干旱氣候與濕潤氣候過渡帶，降雨和蒸發(fā)時空分布不均勻，多年平均水面蒸發(fā)量800～1 000 mm，多年平均降水量800～1 400 mm[18]。研究區(qū)連續(xù)干旱事件較為突出，干旱持續(xù)性較強，如1976-1988年、1985-1988年、1992-2000年均發(fā)生持續(xù)性干旱[19]。

流域雨量站點空間分布如圖2所示。

圖2 息縣流域雨量站點空間分布圖Fig.2 Spatial distribution of rainfall stations in Xixian Basin

3.2 干旱特征變量邊緣分布函數(shù)的推求

淮河上游流域3個干旱特征變量邊緣分布的擬合結果及3個變量間的相關系數(shù)如表1所示。

表1 單變量邊緣分布參數(shù)估計擬合檢驗及變量間相關系數(shù)Tab.1 Fitting trest for estimating parameters of univariate edge distribution and coefficient of correlation among variables

由表1可知，在顯著性水平α=0.05的條件下，淮河上游流域干旱歷時邊緣分布的擬合檢驗中只有三參數(shù)伽馬分布D值小于臨界值，P值大于0.05，能通過K-S檢驗，干旱烈度及干旱烈度峰值的擬合檢驗中只有三參數(shù)韋伯分布能通過K-S檢驗。三個變量之間的相關系數(shù)均大于0.5，表示三個變量間存在較高的兩兩相關性，可構建Copula函數(shù)進行分析。

3.3 多維干旱變量聯(lián)合分布的Copula函數(shù)的推求

3.3.1 二維Copula函數(shù)模型

二維Copula函數(shù)的參數(shù)及擬合優(yōu)度檢驗結果如表2所示。

如表2所示，F(xiàn)rank Copula函數(shù)的擬合優(yōu)度檢驗值RMSE和AIC在4種Copula函數(shù)中最小，該函數(shù)對淮河上游流域二維干旱特征變量聯(lián)合分布的擬合效果最佳。

將最優(yōu)擬合函數(shù)的理論頻率與經驗頻率進行比較，結果如圖3所示。

由圖3可知，以最優(yōu)擬合函數(shù)的理論頻率為縱坐標，對應的經驗頻率為橫坐標所構成的點均勻分布在y=x的兩側，說明擬合效果較好。

表2 二維Copula函數(shù)參數(shù)值及擬合優(yōu)度檢驗Tab.2 Two-dimensional Copula function parameter value and goodness-of-fit test

圖3 二維聯(lián)合分布理論頻率與經驗頻率對比Fig.3 Comparison of theoretical and empirical frequencies of two-dimensional joint distribution

3.3.2 三維Copula函數(shù)模型

Copula函數(shù)在干旱研究中的應用多集中在建立干旱歷時和干旱烈度的二維聯(lián)合分布模型，但干旱特征變量還包括干旱烈度峰值及其他多個相關變量如干旱間隔等[20]，二維Copula函數(shù)模型已不能全面反映干旱的特征。干旱烈度峰值能反映一次干旱事件中最枯月份的干旱程度，因此研究構建了干旱歷時、干旱烈度、干旱烈度峰值三維聯(lián)合分布，結合多個特征變量綜合分析干旱特征，保證研究結果的全面性[2]。

三維Copula函數(shù)的參數(shù)及擬合優(yōu)度檢驗結果如表3所示。

表3 三維Copula函數(shù)參數(shù)值及擬合優(yōu)度檢驗Tab.3 Three-dimensional Copula function parameter value and goodness-of-fit test

由表3可知，F(xiàn)rank Copula函數(shù)的擬合優(yōu)度檢驗值RMSE和AIC在4種Copula函數(shù)中最小，該函數(shù)對淮河上游流域三維干旱特征變量聯(lián)合分布的擬合效果最佳。

將最優(yōu)擬合函數(shù)的理論頻率與經驗頻率進行比較，結果如圖4所示。

圖4 三維聯(lián)合分布理論頻率與經驗頻率對比Fig.4 Comparison of theoretical and empirical frequencies of three-dimensional joint distribution

由圖4可知，以最優(yōu)擬合函數(shù)的理論頻率為縱坐標，對應的經驗頻率為橫坐標所構成的點均勻分布在y=x的兩側，說明擬合效果較好。

3.4 淮河上游干旱重現(xiàn)期計算

依據重現(xiàn)期理論，可以得到邊緣分布為u的單變量重現(xiàn)期計算公式[21]為：

(5)

式中：N是資料的系列長度；m是干旱事件次數(shù)。

多維變量聯(lián)合分布組合重現(xiàn)期包括聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期，邊緣分布分別為u、v的二維聯(lián)合分布聯(lián)合重現(xiàn)期Ta和同現(xiàn)重現(xiàn)期To計算公式[22]為：

(6)

(7)

邊緣分布為u、v和w的三維聯(lián)合分布聯(lián)合重現(xiàn)期Ta和同現(xiàn)重現(xiàn)期To計算公式[22]為：

(8)

To=

(9)

分別計算在單變量重現(xiàn)期為2、5、10、20、50和100 a條件下3個干旱特征變量的邊緣分布函數(shù)值。根據邊緣分布的反函數(shù)推求對應的3個干旱特征變量值。由干旱特征變量值計算最優(yōu)擬合二維及三維Copula函數(shù)值。將Copula函數(shù)值代入式(6)、式(7)以及式(8)、式(9)求出相應的多維組合重現(xiàn)期，并與單變量重現(xiàn)期進行比較，結果如表4所示。

由表4數(shù)據分析可知，隨著單變量重現(xiàn)期的增加，聯(lián)合重現(xiàn)期與同現(xiàn)重現(xiàn)期相應增大且同現(xiàn)重現(xiàn)期的增長幅度遠大于聯(lián)合重現(xiàn)期，在單變量重現(xiàn)期大于20年時尤為明顯。單變量重現(xiàn)期T與多維聯(lián)合分布的聯(lián)合重現(xiàn)期Ta和同現(xiàn)重現(xiàn)期To存在明確的大小關系，從小到大依次為Ta、T、To，組合重現(xiàn)可看作是單變量重現(xiàn)期的兩種極端情況，可作為同水平條件下單變量重現(xiàn)期的閾值，以此可推求干旱特征聯(lián)合概率。在T一定的條件下，所有二維聯(lián)合分布中D與S的Ta最大，對應的To最小，而D與P則相反。這說明D與S的相關性較高，D與P的相關性較低。變量的相關程度越高表明當干旱事件發(fā)生時，兩者同時符合干旱發(fā)生條件的幾率越大。以單變量重現(xiàn)期100年為例，對應的干旱歷時為 14.41個月、干旱烈度為21.91、干旱烈度峰值為3.65。流域發(fā)生干旱歷時大于14.41個月且干旱烈度強于21.91的干旱事件的重現(xiàn)期為482年左右，發(fā)生干旱歷時大于14.41個月且干旱烈度峰值高于3.65的干旱事件的重現(xiàn)期卻達到了1335年左右。流域發(fā)生同時滿足干旱歷時大于14.41個月、干旱烈度強于21.91且干旱烈度峰值高于3.65的干旱事件的重現(xiàn)期已經超過一萬年。說明淮河上游息縣以上流域發(fā)生長歷時強烈度高烈度峰值干旱事件的概率很低。綜上，研究區(qū)長歷時干旱事件往往伴隨著強烈度的特點，但烈度峰值并不會很高。

表4 邊緣分布重現(xiàn)期及多維聯(lián)合分布組合重現(xiàn)期Tab.4 The return period of marginal distribution and the combined return period of multidimensional joint distribution

4 結 語

(1)三參數(shù)伽馬分布對淮河上游息縣以上流域干旱歷時D邊緣分布的擬合效果最好，三參數(shù)韋伯分布對干旱烈度S以及干旱烈度峰值P邊緣分布的擬合效果最好。

(2)Frank Copula函數(shù)對淮河上游息縣以上流域二維及三維干旱特征變量聯(lián)合分布的擬合效果最好。由擬合后函數(shù)的經驗頻率與理論頻率的對比圖可知，函數(shù)點均勻分布在直線 兩側，擬合效果良好。

(3)組合重現(xiàn)期可看作單變量重現(xiàn)期的兩種極端情況，其中聯(lián)合重現(xiàn)期可作為同水平條件下單變量重現(xiàn)期的下閾值，同現(xiàn)重現(xiàn)期可作為同水平條件下單變量重現(xiàn)期的上閾值。結果可為淮河上游息縣以上流域干旱風險的分析預測提供科學依據。

(4)淮河上游息縣以上流域發(fā)生長歷時強烈度高烈度峰值干旱事件的概率很低。流域長歷時干旱事件往往伴隨強烈度的特點，但烈度峰值并不會很高。結論可為流域防旱抗旱工作提供一定的指導。

□