基于多維Copula函數(shù)的淮河上游氣象干旱特征分析

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

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

0 引 言

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

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

1.1 干旱特征變量的識(shí)別與確定

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

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

識(shí)別的具體步驟為：確定SPI閾值R0、R1和R2，分別為-0.5、-1和0。若單月SPI值小于R1，則判定為一次干旱過程[如圖1(a)]，D為1個(gè)月，S為a的陰影面積，P在數(shù)值上等于S；若連續(xù)n月的SPI值均小于R0，則算做一次干旱過程[如圖1(e)]，D為n，S為n月SPI的和，P為SPI的最小值；若兩次干旱過程的時(shí)間間隔為一個(gè)月且該月的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ù)韋伯分布擬合干旱歷時(shí)、干旱烈度以及干旱烈度峰值的邊緣分布。經(jīng)擬合優(yōu)度檢驗(yàn)后，干旱歷時(shí)D采用三參數(shù)伽馬分布擬合，干旱烈度S以及干旱烈度峰值P的邊緣分布采用三參數(shù)韋伯分布擬合。

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

(1)

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

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

(2)

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

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

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

Copula函數(shù)是連接單變量邊緣分布，形成在[0，1]區(qū)間上服從均勻分布的多元函數(shù)[13]。設(shè)F1，F(xiàn)2，…，F(xiàn)n為連續(xù)隨機(jī)變量，則有唯一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ù)結(jié)構(gòu)相對簡單，計(jì)算方便，可以構(gòu)造出多種形式、適應(yīng)性強(qiáng)的多變量聯(lián)合分布函數(shù)，在水文領(lǐng)域的實(shí)際應(yīng)用中較為廣泛。[15]本文采用阿基米德Copula函數(shù)中的四種：Frank、Clayton、Gumbel以及Joe Copula函數(shù)構(gòu)建多維干旱變量的聯(lián)合分布。參數(shù)θ的估計(jì)采用半?yún)?shù)極大似然估計(jì)法[16]。擬合優(yōu)度的檢驗(yàn)選用均方根誤差準(zhǔn)則(RMSE)和赤池信息量準(zhǔn)則(AIC)[17]， 以RMSE和AIC值最小作為擬合最優(yōu)的檢驗(yàn)標(biāo)準(zhǔn)，公式為：

(4)

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

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

2.1 研究區(qū)概況

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

流域雨量站點(diǎn)空間分布如圖2所示。

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

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

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

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

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

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

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

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

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

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

由圖3可知，以最優(yōu)擬合函數(shù)的理論頻率為縱坐標(biāo)，對應(yīng)的經(jīng)驗(yàn)頻率為橫坐標(biāo)所構(gòu)成的點(diǎn)均勻分布在y=x的兩側(cè)，說明擬合效果較好。

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

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

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

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

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

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

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

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

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

由圖4可知，以最優(yōu)擬合函數(shù)的理論頻率為縱坐標(biāo)，對應(yīng)的經(jīng)驗(yàn)頻率為橫坐標(biāo)所構(gòu)成的點(diǎn)均勻分布在y=x的兩側(cè)，說明擬合效果較好。

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

依據(jù)重現(xiàn)期理論，可以得到邊緣分布為u的單變量重現(xiàn)期計(jì)算公式[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計(jì)算公式[22]為：

(6)

(7)

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

(8)

To=

(9)

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

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

表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 結(jié) 語

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

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

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

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

□