基于R-Vine Copula函數的極端降水聯合分布模型及風險識別

曾文穎，徐明慶，宋松柏，吳昊昊

(1.西北農林科技大學水利與建筑工程學院，陜西 楊凌 712100； 2.西北農林科技大學旱區農業水土工程教育部重點實驗室，陜西 楊凌 712100； 3.西北農林科技大學經濟管理學院，陜西 楊凌 712100)

全球氣候變暖加速水文循環，引起降水、蒸發、對流層水汽、徑流等要素發生變化[1]，導致洪水、干旱和熱帶風暴等極端水文事件頻發和加劇[2-3]。由于極端降水突發性強、破壞性大，對社會安全及自然環境造成嚴重影響[4]，尤其對人口密集的城市，暴雨、大雪和冰雹等極端降水天氣影響更大，如2021年河南特大暴雨事件作為強降水致災的典型，對人民生命和財產造成了巨大損失。在氣象監測中有效識別致災性極端降水，對雨洪管理、工程設防標準、災情預警預報等具有重要意義[5]。

風險事件發生的概率常用數理統計法和水文頻率法計算。Smith等[6-7]通過描述連續空間內極值事件的相關性，發現Max-stable模型可以很好地揭示水文氣象極值事件的空間結構與特征；陳星任等[8]采用Pearson相關分析和雙側t檢驗，分析了1961—2018年中國持續極端降水事件的發生頻次、持續日數、總降水量以及極值的變化規律，探討了環流因素對這些變化的影響，發現中國持續極端降水事件的頻率、強度及其對總降水的貢獻率均呈現出增大的趨勢；鄭江禹等[9]從降水量級、降水頻率及發生時間等方面系統分析了珠江流域年、雨季及旱季3個時間尺度上的極端降水特征，探討了極端降水變化特征的機理；劉曾美等[10]采用定性分析法和概率風險分析法，分析了中珠聯圍暴雨與承泄區磨刀門水道的上游西江馬口站的洪水遭遇規律；王輝等[11]對昆明市極端降水事件時空演變特征、嚴重程度及城市效應進行了定量分析，發現昆明市極端降水嚴重程度以及城市化對其影響效應將不斷加劇。極端降水作為水文極值事件，內部屬性要素互相關聯[12]，計算極端降水多因子聯合概率分布，刻畫指標間存在的不同尾部相關性，深入分析相關性規律，有助于提高災害綜合危險性評估準確度。

早期水文的多變量分析主要采用多變量正態分布、對數正態分布、Gumbel分布等方法[13]，這些傳統邊緣分布分析方法通?；谧兞烤€性相關的假設，而實際水文變量可能服從正態或偏態分布，準確構建非線性與非正態的多變量聯合分布是當前研究的熱點和難點[14]。Copula函數可以提供可靠的多變量邊緣分布連接方式，被de Michele等[15]引入降水相關研究，而后Zellou等[16]基于雙變量Joe Copula模型對降水和潮位的遭遇風險及城市危險性進行了分析。但當考慮更多指標及其相關性結構時，Copula函數的參數求解隨建模維度增加變得更加復雜，缺乏靈活性和通用性。

為解決高維Copula函數的求解問題，Bedford等[17]引入Vine圖形化思想，通過構建多元條件相依結構，利用Pair Copula函數將多元聯合分布分解，從而建立Vine Copula模型。常見的C-Vine適合描述有主導變量的相關結構，D-Vine適合描述有臨近關系的結構，均需要對變量之間的相關性作嚴格假設[18]，遍歷性弱[19]。如何建立適合于高維隨機變量的Vine Copula函數來評估極端降水因子之間的聯合分布情況并進行風險分析，找到敏感因子及其對極端降水事件的聯合響應，還需要進一步研究。

本文選用靈活且條件獨立要求弱化的R-Vine Copula函數[20]構建極端降水聯合分布模型進行風險分析，識別致災性敏感因子，并以河南省4個氣象站的實測數據進行模型驗證和對比分析，以期對極端降水事件的風險管控提供參考。

1 模型構建

1.1 邊緣分布函數

在進行極端降水指標邊緣分布擬合時，為消除評價指標之間的量綱差異，需對數據進行歸一化處理[21]。選取常用的伽馬分布(Gamma)、正態分布(Normal)、韋伯分布(Weibull)、廣義極值分布(generalized extreme value，GEV)和對數邏輯分布(log-logistic，LL)5種傳統分布函數[22-23]，以極大似然法(MLE)估計參數，通過擬合優度檢驗，選出描述極端降水指標分布規律最優分布線型，建立邊緣分布。

1.2 R-Vine Copula函數

Copula函數作為邊緣分布函數與多維聯合分布函數之間的紐帶，是描述變量間相依性的一種有效工具[24]。假設u=FX(x)和v=FY(y)分別為連續隨機變量X和Y的累計概率分布函數，由Sklar定理知，隨機變量X和Y的Copula聯合分布函數為[25]

P(X≤x,Y≤y)=C(FX(x),FY(y);θ)=

C(u,v;θ)

(1)

式中θ為衡量邊緣分布之間相關性的Copula參數。

常見的Copula函數類型有Gaussian、Student-t、Clayton、Frank、Gumbel 和Joe等[26]。不同類型Copula函數適用于擬合不同尾部分布特性的序列，其中Gaussian Copula函數尾部漸近獨立，具備零尾部相關性；Student-tCopula函數具備較厚尾部特性，對變量間尾部相關關系較為敏感；Joe Copula函數具有明顯的下尾相關性；Frank Copula函數適合刻畫具有對稱尾部相關性的變量；Gumbel Copula函數對變量分布的上尾部變化比較敏感，能夠捕捉到上尾相關的變化；Clayton Copula函數對變量分布的下尾部變化比較敏感，能夠捕捉到下尾相關的變化。

綜合不同Copula函數的特點，結合圖論和條件函數，建立由節點、樹和邊組成且具備復雜多因子聯合分析能力的R-Vine Copula模型，其中每個節點代表一個變量yi，節點間的連線稱為邊，代表由所連接2個節點組成的Pair Copula函數，后續每一棵樹的節點都來自其前一棵樹的邊[27]，則n維R-Vine Copula函數的密度函數為[28]

(2)

式中：cje,ke|De為邊e對應的Pair Copula密度函數；yDe為由條件集De決定的條件向量；yje、yke分別為節點je和ke處的變量；F(yje|yDe)、F(yke|yDe)為條件分布函數；n為節點數；e為任意邊；Ei為邊的集合；je、ke分別為邊e的兩個節點；De為je和ke間的連接條件集。

當維數較小時可以遍歷搜索所有Vine結構，但隨著維數的升高，Vine結構呈指數增長，求解難度大、可操作性不強[29]。為此，結合TSP(traveling salesperson problem)[30]算法，調用R語言中的TSP包，以每級樹赤池信息準則(AIC)之和最小為依據，優選每個節點處的Copula函數類型并計算參數值。

1.3 基于R-Vine Copula函數的極端降水聯合分布模型

基于R-Vine Copula函數的極端降水聯合分布模型(以下簡稱“R-Vine Copula模型”)構建及風險識別的主要步驟為：①分別建立各極端降水指標的邊緣分布函數，根據KS檢驗、CvM檢驗和均方根誤差進行邊緣分布函數優選，計算原始樣本點處的函數值；②將所選指標兩兩組合成二維數組，以Gaussian Copula、Student-tCopula和Frank Copula，以及Clayton Copula、Gumbel Copula和Joe Copula的0°、90°、180°、270°旋轉式(記Clayton Copula的0°旋轉式為Clayton0°)，共15種Copula函數共同構建一級樹；③以AIC加和值最小為定量評價指標，通過TSP算法和極大似然法選取一級樹的最優Copula函數并計算其參數；④二級及以上的拓撲結構將受到前一級的制約，即上一級樹的一個邊對應于下一級樹的一個節點，重復建立一級樹，直至選出二級樹的最優Pair Copula類型及其參數；⑤重復上述過程，生成所有樹的拓撲結構，完成R-Vine Copula模型的構建；⑥以原始數據均值作為基準值，分別按照不同概率分布(0.1、0.3、0.5、0.7、0.9)，改變其中一個極端降水因子的值，計算聯合概率密度值，進行風險識別。

2 模型驗證和對比分析

2.1 研究區域與數據來源

河南省年平均降水量為600～1 200 mm，淮河以南為1 000～1 200 mm，黃淮之間為700～900 mm，豫北地區為600～700 mm，年降水量區域分布不均，具有從東南向西北減少的趨勢。溫帶大陸性季風氣候的不穩定性，造成年際降水量差別明顯，主要表現為極值比大、變差系數較大以及年際豐枯變化頻繁等。選取河南省鄭州、新鄉、駐馬店和南陽4個氣象站實測數據對模型進行驗證和對比分析。根據中國氣象數據網的中國地面氣候資料日值數據集，搜集4個站點1958—2017年日降水量數據，計算國際通用的4個極端降水指標(大于1 mm年降水總量PT、年降水強度PS(總降水量與降水日數比值)、最大1 d降水量R1和連續5 d最大降水量R5)，繪制各氣象站點PT、R1和R5序列如圖1所示。1958—2017年鄭州站年平均降水量為626.65 mm，最高可達1 028.2 mm；新鄉站在2016年出現日降水量達414 mm的特大降水；駐馬店站年平均降水量為948.78 mm，最高可達1 780.3 mm，1975年和1982年最高日降水量大于400 mm，降水量充沛；南陽站年平均降水量為769.07 mm，最高可達1 434.5 mm?？傮w來看，年降水量在空間分布上由東南向西北逐漸減少，且南部暴雨發生概率高于北部。

(a) 鄭州站

2.2 邊緣分布優選

以水文中常用的Gamma、Normal、Weibull、GEV和LL 5種不同類型分布函數，對歸一化處理后的極端降水指標進行經驗分布和理論分布擬合，采用極大似然法估計參數，通過均方根誤差(RMSE)和KS檢驗進行邊緣分布優選。在優選邊緣分布時，以假設檢驗的概率P>0.05、兩個分布間最大距離D<0.16為臨界值，判斷分布是否可行。以P值接近1、D和RMSE最小為優選標準，綜合考慮3個指標選出不同氣象站不同指標的最優分布函數，結果如表1所示。由表1可見，除鄭州站和駐馬店站的PT和PS以及南陽站的R1指標外，其余均服從GEV分布；4個氣象站R1和R5指標中有87.5%服從GEV分布。由于GEV分布從原數據極值中抽樣，而非依賴原數據概率分布特征，對具有極值的氣象數據邊緣分布擬合較為合理。

表1 極端降水指標邊緣分布函數優選Table 1 Optimization of edge distribution function of extreme precipitation indices

2.3 聯合分布的建立

由優選的邊緣分布，計算樣本點的函數值，代入R-Vine Copula模型中，求解各個連接點處的最優Copula函數類型及其參數，通過AIC值、KS檢驗和CvM檢驗與C-Vine Copula模型(基于C-Vine Copula函數構建的極端降水聯合分布模型)進行對比，結果如表2所示。R-Vine Copula和C-Vine Copula模型的KS檢驗P值均大于0.89，CvM檢驗的P值都接近1，能較好地擬合4個變量之間的聯合分布，模型具有合理性。對比兩種模型的優劣發現，對于鄭州站和駐馬店站，R-Vine Copula模型的KS檢驗更優，南陽站R-Vine Copula模型無論KS檢驗還是AIC值均優于C-Vine Copula模型。整體而言R-Vine Copula模型擬合效果更佳，充分描述了不同變量之間的尾部特征，且無須對數據相依性結構做出預先假設，具有更高的靈活性。

將求解的R-Vine Copula模型各節點Pair Copula函數類型和參數列于表3，以南陽站為例，建立的一級樹中，PT和PS、PT和R1選用Gaussian Copula連接，變量之間漸近獨立，呈現零尾部相關；R1和R5之間具有上尾相關特性，R1的增減直接影響到R5，在一級樹的基礎上，利用最大生成樹原理，繼續生成二級樹，其中“R1,PS;PT”和“R5,PT;R1”之間主要為對稱尾部相關，通過Frank Copula函數建立三級樹。根據計算結果繪制R-Vine結構示意圖如圖2所示。

對R-Vine Copula模型進行檢驗，隨機生成100組數據，計算PT、PS、R1、R5兩兩之間的Kendall和Spearman相關系數，以新鄉站和駐馬店站為例，繪制箱體圖如圖3和圖4所示。新鄉站模擬產生的PT-PS、PT-R1、PS-R1、R1-R5的Kendall相關系數的箱圖中線與原序列接近，其他因子間的Spearman相關系數也分布在箱體的上下四分位數之間。而駐馬店站除PS-R1在箱體上下誤差線之間，其他指標均落在箱體的上下四分位數之間。對圖3進行分析可以發現，無論是新鄉站還是駐馬店站，PT-PS和R1-R5之間具有較大的相關性，可以認為年降水量與降水強度、最大1d降水量與最大5 d降水量的概率分布密切相關，在發生氣象災害時，會出現聯合響應?？傮w而言，R-Vine Copula模型模擬生成的樣本較好地保持了原序列的Kendall相關系數和Spearman相關系數，擬合結果貼近實測值，因此所構建的模型可用于分析研究區極端降水指標的聯合概率分布情況。

表2 R-Vine Copula和C-Vine Copula模型檢驗結果Table 2 Test results of R-Vine Copula and C-Vine Copula model

表3 R-Vine Copula模型各節點函數選擇和參數估計Table 3 Function selection and parameter estimation at each node of R-Vine Copula model

圖2 南陽站R-Vine結構及參數示意圖Fig.2 Schematic diagram of R-Vine structure and parameters at Nanyang station

2.4 風險識別分析

為分析極端降水指標PT、PS、R1、R5對極端降水事件聯合風險的影響，以原始序列的均值作為基準值，采用R-Vine Copula模型計算不同概率分布(0.1、0.3、0.5、0.7、0.9)對應指標值，結果見表4。

在各因子處于同一概率分布時，其對應聯合概率密度值越大，說明該因子對聯合分布變化率的影響越大，即聯合風險對其越敏感[26]。通過改變關鍵風險因子的值，用R-Vine Copula模型計算聯合概率密度的具體方法為：4個極端降水因子中，以其中之一作為關鍵風險因子，另3個因子為非關鍵風險因子，關鍵風險因子取概率分布為0.1、0.3、0.5、0.7和0.9時的統計值，非關鍵風險因子取實測序列的均值來計算聯合概率密度，結果如圖5所示。對于鄭州站，當R5取概率分布0.3、R1取0.5時，對應的聯合概率密度值較大，認為鄭州站的極端降水多因子聯合概率密度主要受到R1和R5的影響，其對聯合分布變化率的影響最大。對于新鄉站和南陽站，PS的概率分布取0.1時，對應的聯合概率密度值最大，概率分布從0.1變化到0.9的過程中，PS從最敏感因子變為敏感性最弱的因子，因此在工程施工安全設計及極端降水災害防范過程中，概率分布為0.1～0.3的PS值應重點考慮。對于駐馬店站，各因子對聯合分布產生影響主要集中在概率分布為0.1～0.3之間，其中PS的概率分布取0.3時，對聯合分布變化率影響最大。

(a) Kendall相關系數

(a) Kendall相關系數

表4 不同概率分布的極端降水指標值Table 4 Extreme precipitation indices with different probability distributions

根據風險識別結果，得出對鄭州站聯合概率密度影響最大的因子為R1和R5，新鄉站、駐馬店站和南陽站則是PT和PS。圖6是其他指標為50%概率分布時，聯合概率密度對敏感因子的聯合響應分布圖，可見聯合概率密度隨敏感因子的概率分布值增大而增大，在敏感因子概率分布值為0.4～0.6附近達到最大，后又隨敏感因子的概率分布值增大而減小，兩個敏感因子對聯合概率密度的影響呈現正相關關系。

3 結 語

本文構建了基于R-Vine Copula函數的極端降水聯合分布模型，該模型無須對數據的相關性結構做出預先假設，具有更高的靈活性。河南省極端降水聯合分布模擬對比結果表明，該模型可以描述變量之間不同的尾部特征，較好地保持了原序列的Kendall和Spearman相關系數等統計特征，可用于分析極端降水指標的聯合概率分布。應用該模型識別了河南省極端降水風險，結果表明：年降水總量PT與年降水強度PS、最大1 d降水量R1與最大5 d降水量R5概率分布呈現密切相關關系；鄭州站的極端降水多因子聯合概率密度主要受R5和R1影響；新鄉站和南陽站的PS隨概率分布增大由最敏感因子降為最不敏感因子；駐馬店站各因子對聯合分布產生影響主要集中在概率分布為0.1～0.3時，且PS的概率分布取0.3時影響最大。

(a) 鄭州站

(a) 鄭州站