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

段奕多 金琳

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

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

中圖分類號 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

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

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

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

1 研究材料與方法

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

1.2 研究方法

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

F（x1，x2，…xn）=C[FX1（x1），F(xiàn)X2（x2），…FXn（xn）] （1）

式中：x1，x2，…xn為觀測樣本，F(xiàn)x為邊緣分布函數(shù)，C為n維Copula函數(shù)

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

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

計算方法：

1.似然值

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

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

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

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

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

3.納什效率系數(shù)（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]為聯(lián)合的經(jīng)驗頻率，[Pi]為聯(lián)合分布的理論頻率，mi為實測樣本中滿足D≤di且S≤si的個數(shù)。

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

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

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

似然函數(shù)：

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）

參數(shù)后驗分布：

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

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

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

2 結(jié)果分析

2.1 Copula函數(shù)

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

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

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

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

2.2 干旱特征分析

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

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

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

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

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

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代表干旱歷時、干旱烈度概率分布函數(shù)

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

3 結(jié)論與討論

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

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

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

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

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

參考文獻

[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]劉章君，郭生練，許新發(fā)，等.Copula函數(shù)在水文水資源中的研究進展與述評[J].水科學進展，2021，32（01）：148-159.

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

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

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

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

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

[10]王鵬新，馮明悅，孫輝濤，等.基于主成分分析和Copula函數(shù)的干旱影響評估研究[J].農(nóng)業(yè)機械學報，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]周平，周玉良，金菊良，等.水文雙變量重現(xiàn)期分析及在干旱中應用[J].水科學進展，2019，30（03）：382-391.

（責編：王慧晴）