基于Copula函數的華南臺風災情重現期研究*

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

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

1 引 言

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

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

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

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

2 理論與方法介紹

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

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

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

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

2.2 Copula函數

2.2.1 Copula函數的定義

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

若函數C為二維Copula函數,在[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為分布函數.

設F是一個三維分布函數，它的邊緣分布可用F1，F2，F3表示.存在一個三維Copula函數C，使得對任意的x1,x2,x3∈R都滿足

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

(1)

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

2.2.2 Copula函數的形式選取及參數估計

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

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

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

2.3 重現期

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

其中，定義單個變量的重現期為

(2)

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

(3)

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

(4)

3 實例分析

3.1 數據來源及因子選取

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

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

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

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

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

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

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

受災人口X1/萬人

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

直接經濟損失X3/億元

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

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

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

農作物受災面積X2/萬公頃

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

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

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

表2 變量間的相關系數表

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

選用Gumbel-Hougaard Copula函數、Clayton Copula函數以及Frank Copula函數來構造兩變量間的聯合分布函數.

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

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

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

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

圖4 受災人口與農作物受災面積的聯合分布函數圖

圖5 農作物受災面積與直接經濟損失的聯合分布函數圖

圖6 受災人口與直接經濟損失的聯合分布函數圖

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

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

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

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

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

3.5 災情重現期

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

表5 三變量在不同重現期下的設計值

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

3.6 災情重現期與臺風災情的比較分析

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

表6 不同重現期范圍的臺風災害個例

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

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

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

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

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

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

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

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

3.7 不同重現期的對比分析

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

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

圖7 不同重現期范圍下的致災臺風與災情臺風

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

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

4 結 論

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

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

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