基于貝葉斯框架下copula函數的干旱特征分析

段奕多 金琳

摘 要：以汾河流域為研究對象，基于SPI采用游處理論提取干旱特征，采用copula函數模擬干旱特征。結果表明：參數值為2.05的Gumbel copula為模擬汾河流域干旱歷時和干旱烈度相依結構的最優copula函數;通過比較擬合效果，MCMC模擬參數估計的精度明顯優于常用的局部優化方法;汾河流域秋冬季干旱頻繁，年際變化呈現旱澇交替趨勢，在研究期（1960—2010年）內共發生54次干旱事件，其中最嚴重的干旱事件重現期為90年一遇，該研究可為汾河流域的干旱預警機制提供參考。

關鍵詞：汾河流域;貝葉斯框架;copula函數;干旱頻率

中圖分類號 P426 文獻標識碼 A文章編號 1007-7731（2021）15-0185-06

Analysis of Drought Characteristics based on Copula in Bayesian Framework

DUAN Yiduo et al.

（College of Water Conservancy and Civil Engineering， Northwest A & F University， Yangling 712100， China）

Abstract： In this paper， Fenhe River Basin is taken as the research object. Based on SPI， drought characteristics are extracted by using the theory of run and stimulated by Copula. The results show that the Gumbel Copula with a parameter value of 2.05 is the best copula function to simulate the correlation between drought duration and drought intensity in Fenhe River Basin; by comparing the fitting results， the accuracy of MCMC simulation parameter estimation is obviously better than the common local optimization method; during the study period， there are total 54 droughts occurred and the most severe drought occurred once in 90 years in Fenhe River Basin. This study can provide a reference for the drought early warning mechanism in Fenhe River Basin.

Key words： Fenhe River Basin; Bayesian framework; Copulas; Drought Frequency

近年來，受全球變暖等影響，氣候和環境發生了重大變化[1]，干旱事件發生頻率增加，造成的危害程度也增大。受地理環境和氣候特征影響，我國是一個干旱災害頻發的國家[2]，據統計，在過去的不到30年中，我國因干旱導致的糧食減產達252億kg，因此對于干旱事件的研究十分重要。在干旱災害的研究內容中，旱災管理一直是學者們重點關注的方面，準確識別干旱事件，并分析干旱事件的重現期，有助于提高農業抗旱的能力，從而減輕農業受干旱災害的影響。

干旱事件是時空連續的具有多維特征的事件，且各個變量之間往往有密切聯系[3]，只針對單變量的統計分析難以提供全面的水文信息，基于此，copula函數被引入對于干旱特征的多元分析中[4]。多年來，copula函數在水文和氣候方面的研究中已經得到了廣泛應用，例如多變量水文事件的頻率分析、水文事件組合分析、水文預報等方面，但大多數水文研究往往局限于固定幾個copula（即Gaussian、t、Frank、Gumbel和Clayton）[4-7]。研究人員從以常用的copula作為備選，進行擬合優度評價，從中選出擬合效果最優的一個，這就使得copula函數的運用在選擇方面受到了局限。在實際中，采用不同的copula函數，計算結果會有很大的差異，選擇最合適的copula函數是運用copula函數分析干旱頻率的基礎，在研究過程中研究者應盡量擴大對模型的研究范圍，嘗試分析新的模型。此外，對于copula函數參數的估計，目前大多數學者采用極大似然法[7，8]，這種方法簡單有效，但容易出現局部最優的情況，并且也無法定量計算出估算頻率的不確定性區間，尤其是在以資料年限短的資料為對象進行分析的情況下，這種方法得出的分析結果可靠性難以判斷。本文采用一種貝葉斯框架下的copula函數方法建立干旱歷時-干旱烈度的聯合概率分布，基于多元copula分析工具箱（Multivariate Copula analysis toolbox，MvCAT），考慮25種copula函數，包括現有水文文獻中沒有充分考慮過的copula函數種類，在貝葉斯框架下使用馬爾科夫鏈蒙特卡洛模估計copula函數的參數，繪制copula的后驗參數分布并估計參數的不確定性，比較各項擬合評價指標，討論得出最優copula函數，從而基于此分析流域多年來的干旱情況。

汾河流域是山西省工業集中、農業發達的主要地區，流經省內6市29縣，流域面積占全省面積的25.5%，水資源總量占全省水資源量的27.2%，汾河是山西省省內第1大河，汾河流域的農業和工業生產值占山西省的50%以上。近年來，汾河流域干旱事件頻發，水文干旱情勢嚴峻[9]，研究汾河流域的干旱特征對于山西省水資源優化配置、流域生態管理等具有一定的參考意義。

1 研究材料與方法

1.1 數據來源 為收集識別干旱事件的樣本，本研究選用標準降雨指數（Standardized Precipitation Index）SPI作為干旱指標，根據游程理論識別干旱事件的歷時和烈度。基礎數據為汾河流域1960—2020年的逐月降雨值，來源于國家氣象科學數據中心（http：//data.cma.cn），資料完整，具有可靠性。選取的水文站為：太原站、介休站、臨汾站和侯馬站，分別位于流域的上、中、下游、中下游。為保證SPI計算準確，資料時間尺度跨越了60年，代表性得到保證，在此時間序列中共識別到54次干旱事件。

1.2 研究方法

1.2.1 Copula函數 Copula是一種構建多個隨機變量聯合分布的函數方法[10]，在[0，1]區間內服從均勻分布。使用Copula函數分析雙變量的聯合分布時，首先要建立單變量概率分布函數，再通過Copula函數把隨機向量X1，X2……Xn的聯合分布函數F（x1，x2，…xn）與各自的邊緣分布函數FX1（x1），…，FXn（xn）相連接的連接函數，即函數C（u1，u2，…，un），使得

F（x1，x2，…xn）=C[FX1（x1），FX2（x2），…FXn（xn）] （1）

式中：x1，x2，…xn為觀測樣本，Fx為邊緣分布函數，C為n維Copula函數

本研究基于多元copula分析工具箱，使用25種copula函數對觀測數據進行擬合并進行擬合優度評價，最終選取排序結果中排名前5的copula函數（Gumbel、Roch-Alegre、BB1、BB5、Tawn）和水文中常用的5個copula（Gaussian、t、Clayton、Frank、Gumbel）進行后續討論。

1.2.2 擬合優度評價 模型擬合優度評價常用的指標有最大似然值、Nash-Sutcliffe效率（NSE）和均方根誤差（RMSE），但這些指標只關注最小化觀測和模型模擬之間的殘差，因此本研究中增加了Akaike信息準則（AIC）和貝葉斯信息準則（BIC），根據式（2）～（8）計算。這些指標中，AIC考慮了模型的復雜性，BIC考慮了模型的復雜度和觀測的數量。RMSE可以量化copula函數模型建立時存在的潛在不確定性，在觀測數據有限的情況下，這對于結果的準確性很重要。模型復雜性越高即模型越靈活，也會更好地擬合觀測數據，但可能會導致模型的過度約束。與似然值相比，AIC考慮了模型的復雜性和誤差殘差的最小化，并提供了1種更穩定的模型預測質量度量，通過根據參數個數增加1個懲罰項，避免了過度約束的問題。

計算方法：

1.似然值

l（θ|[Y]）=[i=1n12πσ2][exp][-12σ-2[yi-yi（θ）]2] （2）

式中：l（θ|[Y]）為似然函數，[Y]為觀測變量的聯合概率分布，σ為誤差標準差的估計值，θ為估計參數，yi為觀測變量

2.貝葉斯信息準則（BIC）

MSE=[1n-1i=1n（Pei-Pi）2] （3）

BIC=nln（MSE）+ln（n） （4）

3.納什效率系數（NSE）

NSE=1-[i=1n（Pei-Pi）2i=1n（Pei-Pi）2] （5）

4.均方根誤差（RMSE）

RMSE=[MSE] （6）

5.信息準則（AIC）

AIC=nln（MSE）+2 （7）

[Pei=mi-0.44n+0.12] （8）

式（2）～（5）中：n為樣本的容量，[Pei]為聯合的經驗頻率，[Pi]為聯合分布的理論頻率，mi為實測樣本中滿足D≤di且S≤si的個數。

AIC、BIC值越低，擬合效果越好。RMSE和NSE的取值范圍分別是RMSE∈[0，∞];NSE=1，NSE∈[-∞，1]，RMSE值越低，擬合效果越好;NSE越接近1，擬合效果越好，函數的參數范圍對它們沒有影響。此外，本文還計算了蒙特卡洛模擬p值，如果p≤0.05，則否定原假設;如果p>0.05，則原假設合理。

1.3 貝葉斯分析 貝葉斯分析是一種模型推理和不確定性量化的有效方法，成功應用于水文學在內的各個領域[11]。當新的觀測樣本被提供時，貝葉斯定理首先計算假設的先驗概率，將模型建立中的所有不確定性歸結為參數，然后根據以下方法估算模型參數的后驗分布。

p（θ|[Y]）=[P（θ）P（Y|θ）P（Y）] （9）

似然函數：

l（θ|[Y]）=[i=1n12πσ2][exp][-12σ-2[yi-yi（θ）]2] （10）

為簡化計算，對式（10）進行簡化為：

l（θ|[Y]）=[--n2ln（2πσ2）-12σ-2i=1n[yi-yi（θ）]2] （11）

參數后驗分布：

P（θ|[Y]）[∝]P（θ）P（[Y]|θ） （12）

式（9）～（12）中：P（θ）和P（θ|[Y]）分別表示參數的先驗分布和后驗分布，σ為觀測誤差的標準查估計值，是水文中最廣泛使用的校準標準[12]，θ為估計參數Y表示copula函數計算的概率值，[Y]表示觀測樣本的聯合概率。P（[Y]|θ）表示似然函數，P（[Y]）為常量，在分析中可以忽略不計。Bayes方程雖然結果準確，但是存在求解困難的局限，因此在計算過程中常使用馬爾科夫鏈蒙特卡羅算法進行后驗分布的取樣。

1.4 馬爾科夫鏈蒙特卡羅模擬 馬爾科夫鏈蒙特卡羅算法（MCMC）是1種從高維的、復雜的分布中進行抽樣的統計方法，用來描述基于貝葉斯框架的后驗參數分布[13]。混合演化MCMC的優勢在于智能起點的選擇，采用自適應采樣算法（Adaptive metropolis，AM）、差分演化算法（Differential evolution，DE）和斯諾克更新算法來計算可行域，實現在1次計算中，同時采用多個起點，搜索多個計算空間，形成多條計算鏈，找到全局最優解的估計值和參數的后驗分布。在計算時，將樣本隨機分配給n個Copula函數，并選擇似然值最高的樣本作為馬爾可夫鏈的起點。與目前copula函數參數估計常用的局部優化算法相比，MCMC算法計算范圍更廣，從多個起點出發因此可以確保估計的參數為全局最優，并可以描述參數的不確定性[13]。

2 結果分析

2.1 Copula函數

2.1.1 擬合優度評價 模型的選擇一般基于擬合優度評價，常用的指標有最大似然值、Nash-Sutcliffe效率（NSE）和均方根誤差（RMSE），除此之外，本研究還選取了Akaike信息準則（AIC）、貝葉斯信息準則（BIC）作為參考，計算結果見表1。

由于不同的指標關注重點不同，根據不同指標選擇的最優模型也不同，研究者應該根據基于copula模型解決問題的目的來選擇。本研究主要關注干旱極端事件，對于模型尾部的擬合效果要求更高，因此確定Gumbel copula來模擬干旱變量之間的關系，聯合概率分布見圖1，函數表達式為：

C（u，v）=exp[-[（-lnu）1θ+（-lnV）1θ]θ] （13）

2.1.2 參數估計及不確定性分析 使用貝葉斯分析和MCMC模擬來尋找參數的后驗分布，為了量化copula模型的不確定性，繪制了一組Copula族的后驗參數分布圖，見圖2。圖中上方的符號表示局部優化方法得到的copula函數的參數，而柱狀圖表示根據MCMC模擬得到的參數，橫坐標上的“×”代表MCMC的最大似然參數。由圖2可以看出，圖中上方的符號與下方的“×”位置明顯不同，表明這2種方法估算的參數有明顯差異。MCMC模擬模型參數的RMSE為0.0763，而局部優化算法模擬參數的RMSE值為0.2231，在擬合度上全局優化算法相對局部優化算法有明顯的優勢。由圖2（b）發現，Gaussian copula的參數分布接近均勻分布，這表明觀測數據包含的信息不足以約束copula的參數。另外，由圖2（c）發現，Tawn copula的第2個參數分布在參數的邊界上得到合并，最大似然參數位于邊界上，這意味著優化算法試圖通過越界來提高擬合度，這是不合理的，也表明選擇這樣的copula對于數據的擬合不恰當。

2.2 干旱特征分析

2.2.1 干旱事件 為了比較不同時間尺度下的干旱特征，本文分別計算了時間尺度為3個月、6個月和12個月的SPI，不同時間尺度的SPI序列變化趨勢基本相同，都反映了流域旱澇交替的特征。時間尺度越小，SPI隨時間的波動頻率越高，時間尺度越大，SPI序列越平坦。本研究選擇3個月尺度的SPI，以反映農業干旱情況[9]。將干旱劃分為特旱、重旱、中旱、輕旱幾個等級，在流域的干旱情況中，秋冬季節在1年中干旱較為嚴重。在研究區內，發生了3次特旱事件，分別在1965年冬季、1965年秋季和1997年秋季;春季只在2000年發生了1次嚴重干旱，夏季在1968年和2001年和發生了嚴重干旱，其余春夏季節發生的干旱事件均集中在輕旱、中旱;秋季分別在1965年和1997年發生了特旱，在1972年和1986年發生了嚴重干旱;冬季在1965年發生了特旱，在1970年、1997和1998年發生了嚴重干旱。在研究的時間序列中，秋冬季發生干旱的頻率也明顯高于春夏。

2.2.2 干旱重現期 干旱事件的頻率往往用重現期來表示，是水文信息的一個重要參數[15]。計算重現期，首先要建立合適的多變量的聯合概率分布，然后才可以根據式（14）-（16）計算。單變量的重現期計算公式為：

T=[L1-F（x）] （14）

式中：T為重現期，F（x）為累積概率分布函數，[L]為干旱事件的平均時間間隔。

對于雙變量干旱事件，在隨機變量的取值范圍內，各個變量同時大于某一設定的值，重現期的計算公式為：

P（T≥t∩S≥s）=1-u-v+C（u，v） （15）

Tds=[LP（T≥t∩S≥s）]=[L1-u-v+C（u，v）] （16）

式中：T、S表示干旱時間、干旱烈度，u、v代表干旱歷時、干旱烈度概率分布函數

用Gumbel copula來模擬干旱歷時和干旱變量之間的依賴結構，分析干旱重現期，干旱變量的重現期見圖3。首先計算干旱間隔時間為13.38（月），根據式（13）～（15）計算得到，對于單變量重現期而言，研究序列中發生的最嚴重干旱事件的重現期（干旱歷時、干旱烈度）分別為9.13年一遇、43.2年一遇。就雙變量聯合重現期而言，在所研究的時間序列中，發生最嚴重的干旱事件，同現重現期為32年一遇、而聯合重現期為90年1遇。比較分析近60年里汾河上中游地區發生的幾次典型干旱事件，1965年發生29年一遇的干旱事件持續了7個月，直到1966年3月旱期才開始減輕。1972年發生干旱同樣連續持續了7個月，干旱程度主要為中旱，這次干旱事件嚴重程度相比1965年程度較輕，但旱情反復出現，在接近1年之后才結束，重現期約為26年一遇。此外，在1986年發生了持續7個月的干旱，干旱程度也均為重旱、特旱;1997年發生干旱持續了6個月的干旱事件，干旱程度幾乎持續保持特旱，重現期為28年一遇。

3 結論與討論

干旱事件的各個變量之間是密切相關的，研究它們之間的相互關系對于干旱預測和風險評估至關重要，copula函數可以很好地模擬干旱特征變量之間的關系、模擬干旱事件的發生，是干旱分析的有效工具，對于旱情的監測和預報有重要意義。本文基于汾河流域近60年的逐月降水數據，以SPI3為指標基于游程理論識別干旱事件，采用多元copula分析工具箱在貝葉斯框架下使用MCMC模擬來估計copula參數及其潛在不確定性，計算擬合copula概率的等值線，再利用Copula函數建立的聯合分布分析汾河流域的干旱特征，得到的結論有：

（1）由SPI很好的反映汾河流域的多年干旱情況：在長時間序列下流域表現出旱澇交替的特點;秋冬季節的干旱明顯比春夏季節嚴重，干季的干旱比濕季嚴重;且觀察多年SPI的變化趨勢，發現整體呈現輕微下降趨勢，證明汾河流域的干旱情勢并不樂觀。

（2）通過干旱識別、擬合優度評價選擇copula函數、估計函數參數等步驟，擬合汾河流域干旱事件歷時和烈度的最優copula為參數值等于2.0461的Gumbel copula，通過干旱頻率分析，汾河流域在1960—2020年內發生最嚴重的干旱為90年一遇;

（3）貝葉斯方法是模型參數估計的有效方法，相對于以往常用的局部優化方法具有一定的優勢。在以往研究中常用的局部優化算法會存在陷入局部最優的現象，導致結果有偏差，而貝葉斯框架下的MCMC模擬，可以得到全局最優解，還能提供參數估計的潛在不確定性。尤其是數據信息有限時，copula函數參數估計存在很大的不確定性。

本研究對于干旱事件的描述僅選取了干旱歷時和干旱烈度2個變量，選擇的特征變量越多，在計算過程中建立多變量間的copula函數越復雜，參數估計也越難計算。因此我們還需要考慮如何提出更簡單有效的方法，在基于多個特征變量的基礎上降低計算的復雜性。此外，目前水文文獻中常用的copula函數僅局限于固定的幾個（Gaussian、t、Frank、Gumbel和Clayton），但存在的問題是這些copula函數有時不能很好的描述變量間的依賴結構。

參考文獻

[1]Zhang Q， Yao Y， Li Y， et al. Causes and Changes of Drought in China： Research Progress and Prospects[J].Journal of Meteorological Research， 2020，34（3）：460-481.

[2]屈艷萍，呂娟，蘇志誠，等.抗旱減災研究綜述及展望[J].水利學報，2018，49（01）：115-125.

[3]徐翔宇，許凱，楊大文，等.多變量干旱事件識別與頻率計算方法[J].水科學進展，2019，30（03）：373-381.

[4]劉章君，郭生練，許新發，等.Copula函數在水文水資源中的研究進展與述評[J].水科學進展，2021，32（01）：148-159.

[5]陳文華，徐娟，李雙成.怒江流域下游地區氣象與水文干旱特征研究[J].北京大學學報（自然科學版），2019，55（04）：764-772.

[6]韓會明，劉喆玥，劉成林，等.基于Copula函數的贛江流域氣象干旱特征分析[J].水電能源科學，2020，38（08）：9-13.

[7]龍瑞昊，暢建霞，張鴻雪，等.基于Copula的瀾滄江流域氣象干旱風險分析[J].北京師范大學學報（自然科學版），2020，56（02）：265-274.

[8]覃麗華，張鑫.基于連續干旱游程的干旱特征分析與應用[J].中國農村水利水電，2020（10）：164-169.

[9]趙雪花，趙茹欣.水文干旱指數在汾河上游的適用性分析[J].水科學進展，2016，27（04）：512-519.

[10]王鵬新，馮明悅，孫輝濤，等.基于主成分分析和Copula函數的干旱影響評估研究[J].農業機械學報，2016，47（09）：334-340.

[11]Prakash A，Panchang V，Ding Y，et al. Sign Constrained Bayesian Inference for Nonstationary Models of Extreme Events[J]. Journal of Waterway Port Coastal and Ocean Engineering， 2020， 146（5）：04020029.

[12]Gupta V，Jain M K，Singh V P. Multivariate Modeling of Projected Drought Frequency and Hazard over India[J]. Journal of Hydrologic Engineering，2020，25（4）：04020003-1-19.

[13]Sadegh M，Ragno E，Aghakouchak A. Multivariate Copula Analysis Toolbox （MvCAT）：Describing dependence and underlying uncertainty using a Bayesian framework[J]. Water Resources Research，2017，53.

[14]Zhu S，Luo X，Chen S，et al. Improved Hidden Markov Model Incorporated with Copula for Probabilistic Seasonal Drought Forecasting[J]. Journal of Hydrologic Engineering，2020，25（6）：04020019.

[15]周平，周玉良，金菊良，等.水文雙變量重現期分析及在干旱中應用[J].水科學進展，2019，30（03）：382-391.

（責編：王慧晴）