基于Copula函數(shù)的華南臺風災情重現(xiàn)期研究*

王 萌，劉合香，盧耀健，李廣桃

(南寧師范大學 數(shù)學與統(tǒng)計科學學院，廣西 南寧 530001)

1 引 言

華南地區(qū)是中國遭受臺風影響最嚴重、最頻繁的地區(qū)之一[1].如2015年臺風“彩虹”、2016年臺風“莎莉嘉”和2017年臺風“天鴿”的登陸，都對華南地區(qū)造成了嚴重的災情影響.研究登陸華南地區(qū)的臺風災害有著比較重要的意義.

近年來，一些學者通過建立不同的模型對華南地區(qū)臺風災害展開了災情研究.陳仕鴻(2011)等[2]根據(jù)廣東省臺風災害歷史數(shù)據(jù)，將選取的5個評估因子分為5個等級，建立了離散型Hopfield神經(jīng)網(wǎng)絡模型，對廣東省的臺風災情進行評估.陳燕璇(2016)等[3]采用主成分分析(PCA)、等距特征映射(ISOMAP)和信息熵特征提取方法，將選取的致災因子和承災體因子作為輸入，災情等級作為輸出，建立了登陸廣西臺風的災情概率神經(jīng)網(wǎng)絡模型，對臺風災情進行預評估.

Copula函數(shù)常被應用于水文災害研究中，為了更準確地研究自然災害的重現(xiàn)期，一些學者通過建立Copula函數(shù)重現(xiàn)期模型來對自然災害的災情進行分析.萬永靜(2017)等[4]采用Copula函數(shù)，構建了不同重現(xiàn)期下的降雨量與潮位的二維聯(lián)合分布函數(shù)，計算兩者的遭遇頻率.馬建琴(2017)等[5]對提取的水文干旱特征變量，利用多種Copula函數(shù)模型擬合其聯(lián)合分布函數(shù)，計算對應的聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期，對水文干旱的發(fā)生頻率進行了分析.

借助Clayton Copula函數(shù)，構造災情因子受災人口、農(nóng)作物受災面積和直接經(jīng)濟損失的三維聯(lián)合分布，計算華南臺風災害的災情重現(xiàn)期，探討臺風災害的發(fā)生頻率，將會是有價值的研究工作.

2 理論與方法介紹

2.1 單變量邊緣分布函數(shù)的選取及參數(shù)估計

利用Copula函數(shù)對選取的變量構造聯(lián)合分布之前，需要確定各變量的邊緣分布函數(shù).選取一些常用的分布來構造邊緣分布.這些分布分別為指數(shù)分布(Exp(λ))、伽馬分布(Ga(α,λ))、正態(tài)分布(N(μ,σ2))、對數(shù)正態(tài)分布(LN(μ,σ2))和韋布爾分布(Weibull(a,b,c)).

使用矩法估計對單變量的不同邊緣分布函數(shù)進行參數(shù)估計.采用矩法估計參數(shù)是因為它具有易于使用和不需要知道總體分布的優(yōu)勢.

選用Kolmogorov-Smirnov檢驗法(簡稱K-S)，對單變量的邊緣分布函數(shù)進行擬合優(yōu)度檢驗，選出最優(yōu)的擬合分布函數(shù).

2.2 Copula函數(shù)

2.2.1 Copula函數(shù)的定義

由Nelsen(2006)[6]可得到二維Copula函數(shù)的定義.

若函數(shù)C為二維Copula函數(shù),在[0,1]2→[0,1]則滿足：

(1)?u,v∈[0,1],有C(u,0)=C(0,v)=0,C(u,1)=C(1,v)=v，

(2)?[u1,u2]×[v1,v2]?[0,1]2,有C(u2,v2)-C(u2,v1)-C(u1,v2)+C(u1,v1)≥0.

其中u,v，u1,u2,v1,v2為分布函數(shù).

設F是一個三維分布函數(shù)，它的邊緣分布可用F1，F(xiàn)2，F(xiàn)3表示.存在一個三維Copula函數(shù)C，使得對任意的x1,x2,x3∈R都滿足

F(x1,x2,x3)=C(F1(x1),F2(x2),F3(x3)).

(1)

其中C:[0,1]3→[0,1]，F(xiàn)1(x1),F2(x2),F3(x3)分別表示變量X1、X2和X3的分布函數(shù).特別地，若F1，F(xiàn)2，F(xiàn)3是連續(xù)函數(shù)，則C是唯一確定的.

2.2.2 Copula函數(shù)的形式選取及參數(shù)估計

在Archimedean(阿基米德)型Copula函數(shù)中，Gumbel-Hougaard(簡稱GH)、Clayton和Frank Copula函數(shù)被廣泛應用于水文以及自然災害研究中[7].選取這三種Copula函數(shù)構建變量的聯(lián)合分布函數(shù).

對二維Copula函數(shù)進行參數(shù)估計時，采用的是相關性指標法[8].該方法需要先計算兩變量的Kendall秩相關系數(shù)τ，再根據(jù)τ與Archimedean Copula函數(shù)的參數(shù)θ之間的關系式求解得到參數(shù)θ的值.

在估計三維Copula函數(shù)的參數(shù)時，選用極大似然法.它不僅克服了運算復雜的問題，而且還能取得較為優(yōu)異的估計結果[9-10].根據(jù)式(1)求出其對應的三變量聯(lián)合分布密度函數(shù)，構造三變量的似然函數(shù)，得到對數(shù)似然函數(shù)，計算對數(shù)似然函數(shù)的最大值點，即可得到三維聯(lián)合分布函數(shù)中參數(shù)θ的最大似然估計值.

2.3 重現(xiàn)期

對于要計算的重現(xiàn)期是指研究災情類似的臺風災害平均多少年重復出現(xiàn)一次.通過選取的臺風災情因子計算得到的重現(xiàn)期來衡量某類臺風災害發(fā)生的頻率.其重現(xiàn)期可用T表示，常用的單位是 “年”.

其中，定義單個變量的重現(xiàn)期為

(2)

三變量聯(lián)合重現(xiàn)期是指變量X1,X2或X3超過某個特定值的重現(xiàn)期，其計算公式為式(3).

(3)

三變量同現(xiàn)重現(xiàn)期是指變量X1,X2和X3都超過某個特定值的重現(xiàn)期，見式(4).

(4)

3 實例分析

3.1 數(shù)據(jù)來源及因子選取

為了描述臺風災害的災情情況，從人口、建筑、農(nóng)業(yè)和經(jīng)濟四個方面選取受災人口、倒塌房屋數(shù)、農(nóng)作物受災面積和直接經(jīng)濟損失作為災情因子進行分析[11].再從臺風災害數(shù)據(jù)的完整性考慮,選取的臺風災情數(shù)據(jù)來源于1981-1999年中國氣象災害大典[12-14](廣東卷、海南卷、廣西卷)、2000-2003年廣東省、海南省和廣西氣候中心和2004-2014年中國熱帶氣旋年鑒[15].2015-2016年災情數(shù)據(jù)來源于廣西氣候中心、廣東氣候中心和海南氣候中心，2017年災情數(shù)據(jù)來源于中國氣象災害年鑒[16].

為了消除各災情因子間的量綱影響，對其原始數(shù)據(jù)進行歸一化處理.利用Spearman和Kendall秩相關系數(shù)法對災情因子進行兩兩間檢驗.得出受災人口、農(nóng)作物受災面積和直接經(jīng)濟損失兩兩因子間具有較強的相關性.從臺風災害的災情因子間的強相關性及代表性考慮，選取臺風災害的受災人口(X1)、農(nóng)作物受災面積(X2)和直接經(jīng)濟損失(X3)三個災情因子來研究華南臺風災害.

3.2 單變量分布的最優(yōu)選取及參數(shù)估計

針對選取的受災人口(X1)、農(nóng)作物受災面積(X2)和直接經(jīng)濟損失(X3)，從常用的指數(shù)分布、伽馬分布、正態(tài)分布、對數(shù)正態(tài)分布和韋布爾分布中，利用矩法估計變量的參數(shù)，結合K-S(Kolmogorov-Smirnov)檢驗法，分別選取這三個變量的最優(yōu)邊緣分布，其中三變量的具體分布情況見表1.

根據(jù)K-S檢驗的原則，其值越小說明該分布擬合的效果越好.從表1中可看出，受災人口的K-S檢驗值最小為0.1918，其對應的分布為正態(tài)分布，說明受災人口用正態(tài)分布擬合最優(yōu).農(nóng)作物受災面積在對數(shù)正態(tài)分布下的檢驗值最小，為0.1422，應該用對數(shù)正態(tài)分布擬合該變量的分布.直接經(jīng)濟損失擬合最優(yōu)的分布為對數(shù)正態(tài)分布，經(jīng)K-S值計算，該變量在對數(shù)正態(tài)分布下的檢驗值為0.1125，是最小值.

表1 不同分布的參數(shù)估計值及擬合優(yōu)度檢驗

根據(jù)X1、X2和X3三變量的擬合分布函數(shù)，可繪出它們各自所服從的最優(yōu)分布函數(shù)與對應的經(jīng)驗分布函數(shù)擬合曲線圖進行對比分析，具體如圖1-3所示.

受災人口X1/萬人

從圖1-圖3中可以看出，隨著單變量觀測值的增加，對應的經(jīng)驗分布函數(shù)值逐漸接近于1.選取的最優(yōu)分布的函數(shù)值也隨著對應的單變量觀測值的增大而越接近于1.

直接經(jīng)濟損失X3/億元

從這3個變量的分布擬合圖中還可以發(fā)現(xiàn)，這三個變量的經(jīng)驗分布函數(shù)曲線圖與其對應的最優(yōu)分布函數(shù)曲線圖相比，兩者的趨勢基本一致，整體都大致呈逐步遞增趨勢.

從圖1可看出受災人口用正態(tài)分布擬合其經(jīng)驗分布，雖然擬合值與實際觀測值相比偏小，但整體趨勢一致，且誤差較小，表明用正態(tài)分布擬合受災人口這一變量較為合適.

從圖2可以看出用對數(shù)正態(tài)分布擬合農(nóng)作物受災面積，與實際觀測值的經(jīng)驗分布函數(shù)曲線大致重合，只有個別觀測值存在誤差，且誤差值較小.根據(jù)兩者分布函數(shù)曲線圖的對比，可得農(nóng)作物受災面積這一變量用對數(shù)正態(tài)分布擬合較為合適.

農(nóng)作物受災面積X2/萬公頃

由圖3可得，直接經(jīng)濟損失的經(jīng)驗分布函數(shù)圖與對數(shù)正態(tài)分布函數(shù)圖的曲線幾乎完全重合，說明可以用對數(shù)正態(tài)分布擬合其分布.

3.3 二維Copula函數(shù)的參數(shù)估計及擬合優(yōu)度檢驗

利用Spearman和Kendall秩相關性系數(shù)檢驗方法，對受災人口(X1)、農(nóng)作物受災面積(X2)和直接經(jīng)濟損失(X3)兩兩變量間進行了相關性系數(shù)檢驗，具體結果見表2.

表2 變量間的相關系數(shù)表

由表2的計算結果可看出，不管是Spearman檢驗法還是Kendall檢驗法，這三個變量兩兩之間的相關系數(shù)都比較大，表明兩兩變量間均具有較強的正相關性，也表明了選取的這三個變量間可通過Copula函數(shù)構造聯(lián)合分布.

選用Gumbel-Hougaard Copula函數(shù)、Clayton Copula函數(shù)以及Frank Copula函數(shù)來構造兩變量間的聯(lián)合分布函數(shù).

對受災人口(X1)、農(nóng)作物受災面積(X2)和直接經(jīng)濟損失(X3)兩兩變量之間構造Copula聯(lián)合分布函數(shù)，利用RMSE和AIC信息準則檢驗法，根據(jù)其值越小擬合效果越好的檢驗準則，選出兩兩變量之間的最優(yōu)聯(lián)合分布函數(shù).其中X1-X2、X2-X3和X1-X3的參數(shù)擬合優(yōu)度評價結果見表3.

從表3可知，變量X1-X2、X2-X3和X1-X3均在Clayton Copula函數(shù)下的RMSE和AIC值最小，其RMSE最小值分別為0.2893、0.3003和0.2953，AIC最小值分別為0.3307、0.2734和-0.1223.結合RMSE和AIC的判斷準則，可得出受災人口與農(nóng)作物受災面積、農(nóng)作物受災面積與直接經(jīng)濟損失、受災人口與直接經(jīng)濟損失變量間都用Clayton Copula函數(shù)擬合效果較好.

表3 二維Copula函數(shù)的參數(shù)及擬合優(yōu)度

根據(jù)X1-X2、X2-X3和X1-X3變量間的最優(yōu)擬合分布函數(shù)，可繪出它們的聯(lián)合分布函數(shù)圖，具體見圖4-圖6.

圖4 受災人口與農(nóng)作物受災面積的聯(lián)合分布函數(shù)圖

圖5 農(nóng)作物受災面積與直接經(jīng)濟損失的聯(lián)合分布函數(shù)圖

圖6 受災人口與直接經(jīng)濟損失的聯(lián)合分布函數(shù)圖

根據(jù)圖4、圖5和圖6可分別得出受災人口和農(nóng)作物受災面積、農(nóng)作物受災面積和直接經(jīng)濟損失及受災人口和直接經(jīng)濟損失的聯(lián)合分布函數(shù)圖.從圖4、圖5和圖6的圖像漸變過程可以看出，顏色越淺，聯(lián)合分布函數(shù)所占的面積越小，聯(lián)合分布的函數(shù)值越接近于1.顏色越深，聯(lián)合分布函數(shù)所占的面積越大，聯(lián)合分布的函數(shù)值越趨于0.如圖4所示，隨著受災人口和農(nóng)作物受災面積變量的分布函數(shù)值增大，其對應的Clayton Copula函數(shù)聯(lián)合分布值就越趨于1.

3.4 三維Copula函數(shù)的參數(shù)估計及擬合優(yōu)度檢驗

采用極大似然法對三維Copula函數(shù)的參數(shù)進行估計，并利用RMSE和AIC信息準則對所計算的參數(shù)進行擬合優(yōu)度檢驗.可得三變量Copula函數(shù)的參數(shù)估計和和擬合優(yōu)度結果，如表4所示.

表4 三維Copula函數(shù)的參數(shù)及擬合優(yōu)度

由以上兩種擬合優(yōu)度檢驗方法的評價準則和表4可得，三變量用Clayton Copula函數(shù)擬合的RMSE值為0.2738、AIC值為-0.0162，這兩個評價準則得出的值都是最小的.表明受災人口、農(nóng)作物受災面積和直接經(jīng)濟損失用Clayton Copula函數(shù)擬合的效果較好.

3.5 災情重現(xiàn)期

利用公式(2)-式(4)，結合式(1)，分別求出不同單變量重現(xiàn)期、聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期下的受災人口(X1)、農(nóng)作物受災面積(X2)和直接經(jīng)濟損失(X3)的設計值，其計算結果見表5.

表5 三變量在不同重現(xiàn)期下的設計值

從表5可以發(fā)現(xiàn)，在聯(lián)合重現(xiàn)期下，受災人口、農(nóng)作物受災面積和直接經(jīng)濟損失的設計值都大于單變量重現(xiàn)期和同現(xiàn)重現(xiàn)期下的設計值，表明在預測臺風災害發(fā)生的頻率時用三變量聯(lián)合重現(xiàn)期預測臺風災害的發(fā)生頻率會更加安全，并將其三變量聯(lián)合重現(xiàn)期記為災情重現(xiàn)期.

3.6 災情重現(xiàn)期與臺風災情的比較分析

表6所示的是不同重現(xiàn)期范圍內(nèi)的49個臺風個案.因災情重現(xiàn)期描述的是某類臺風災害發(fā)生的頻率，通常頻率用P表示，以百分比(%)作單位，頻率與重現(xiàn)期之間的關系可以表示為P=1/T.對于某種臺風災害來說，其重現(xiàn)期越長意味著這類臺風災害發(fā)生的頻率越小.

表6 不同重現(xiàn)期范圍的臺風災害個例

從表6中可以看到，重現(xiàn)期區(qū)間為0～2年時，登陸的臺風最多，有35個臺風.在重現(xiàn)期區(qū)間為2～5年時，有3個臺風登陸.在重現(xiàn)期區(qū)間為5～10年時，登陸的臺風個數(shù)為8個.在重現(xiàn)期區(qū)間為10～20年時，發(fā)生的臺風災害有2個.在重現(xiàn)期區(qū)間為20～50年時，發(fā)生的臺風災害有1個.

說明重現(xiàn)期越短的臺風災害發(fā)生的頻率越大，重現(xiàn)期越長的臺風災害發(fā)生的頻率就越小.

綜合單個臺風的災情重現(xiàn)期可以預測臺風災害發(fā)生的頻率.由表6可知，0518號臺風“達維”的災情因子重現(xiàn)期落在5～10年內(nèi)，那么可預測在2005年后的5～10年內(nèi)還會至少出現(xiàn)一次和0518號臺風“達維”所造成災情影響類似的臺風災害.

根據(jù)聯(lián)合重現(xiàn)期的計算公式可得出，0518號臺風“達維”的重現(xiàn)期為5.2376年，那么可預測在2011年左右會出現(xiàn)一次與0518號臺風“達維”災情類似的臺風災害.

由表6可知，2011年發(fā)生的1117號臺風“納沙”與臺風“達維”造成的災情情況類似.0518號“達維”和1117號“納沙”臺風都為強臺風.

以災情重現(xiàn)期作為判斷依據(jù)，對2015-2017年華南地區(qū)的臺風災情數(shù)據(jù)進行重現(xiàn)期范圍判斷，然后再利用公式(3)得出這三年的臺風的重現(xiàn)期值，最后臺風災害的重現(xiàn)期區(qū)間范圍和重現(xiàn)期值進行對比分析，結果見表7.

表7 2015-2017年華南地區(qū)臺風災害的災情重現(xiàn)期

根據(jù)1117號臺風“納沙”的重現(xiàn)期結果，預測在2017年發(fā)生的臺風災害所造成的災情與其類似.在2017年發(fā)生的1713號臺風“天鴿”所造成的災情影響確實如此.

3.7 不同重現(xiàn)期的對比分析

將表6的臺風災情重現(xiàn)期范圍下的49個臺風與致災重現(xiàn)期[17]范圍下得到的臺風進行對比分析.發(fā)現(xiàn)有個別的臺風災害的致災因子重現(xiàn)期和災情因子重現(xiàn)期落在不同的重現(xiàn)期范圍.造成這一現(xiàn)象的主要原因是因為臺風災害造成的災情不僅受致災因子的影響，還會受到承災體因子的影響、當?shù)仡A報、預防臺風的及時性以及災后救援等客觀因素的影響.

對落在致災重現(xiàn)期的不同重現(xiàn)期范圍內(nèi)的臺風數(shù)目與落在不同災情因子重現(xiàn)期范圍內(nèi)的臺風數(shù)目進行對比分析，具體結果見圖7.

圖7 不同重現(xiàn)期范圍下的致災臺風與災情臺風

從圖7觀察到，在致災重現(xiàn)期與災情重現(xiàn)期的各階段重現(xiàn)期范圍內(nèi)，其對應的臺風個例總數(shù)目相差不大.

綜上所述，與之前致災重現(xiàn)期計算的臺風災害的重現(xiàn)期范圍相比，由于其他的致災因素和客觀因素導致存在個別同一編號的臺風災害的重現(xiàn)期范圍不一致，但整體上大多數(shù)的臺風災害的重現(xiàn)期范圍還是一致的.表明利用臺風災情計算的重現(xiàn)期具有一定的可靠性.為了更好地對以后發(fā)生的臺風災害進行預測，可結合災情重現(xiàn)期和致災重現(xiàn)期進行綜合性分析.

4 結 論

選取華南臺風災害的受災人口、農(nóng)作物受災面積和直接經(jīng)濟損失三個災情影響因子，利用Copula函數(shù)的相關理論，構造三維Clayton Copula聯(lián)合分布函數(shù)，計算災情重現(xiàn)期，分析華南臺風災情，得出以下結論.

對比單變量重現(xiàn)期、聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期下的三變量設計值發(fā)現(xiàn)，采用三變量聯(lián)合重現(xiàn)期來分析臺風災情會更加有效.

通過計算華南地區(qū)發(fā)生的49個臺風的災情重現(xiàn)期，得出臺風災害在重現(xiàn)期0～2年內(nèi)發(fā)生的頻率最高，在重現(xiàn)期為20～50年內(nèi)發(fā)生的頻率最低.