基于Copula函數的降雨隨機模擬方法在古黃河流域的應用

潘奇鑫，李 帆，劉 萍，李春秀

(揚州大學水利科學與工程學院，江蘇 揚州 225009)

在極端氣候頻繁出現的背景下，中國洪澇、干旱災害頻發，防范洪澇、干旱災害的任務變得越來越嚴峻。自唐以來，江蘇地區就是全國經濟發達的地區之一，也是自然災害頻發的地區，尤以洪澇、干旱災害最為顯著[1]。中長期水文預報作為預防災害的有效手段之一，在近些年得到了廣泛的應用，其對水庫合理調度、防洪減災、水資源利用有著重要意義。而現階段的水文預報大多具有確定性，回避了預報值的不確定性，無法滿足對風險信息的需求。以概率的形式定量分析水文預報不確定性，據此做出的概率水文預報，不僅在理論上更加科學，而且在實踐應用中能產生更高的效益，是水文預報技術發展的必然趨勢。張銘等[2]采用貝葉斯概率水文預報系統在中長期徑流預報中提高水文預報的精度。中長期概率預報不僅可得到比確定性預報更高的精度，而且還能為決策者提供更豐富的預報信息[3]。

相對于傳統的水文頻率研究方法，Copula函數法具有“連接”的功能[4]，且計算過程簡便靈活，同時適用于邊緣分布相同與不同的聯合分布，也適用于多維的聯合分布，且效果相對較好，基本能解決傳統水文頻率研究方法中存在的問題。起初，Copula函數主要在金融、保險等領域運用，自2003年起Copula函數才開始在水文頻率計算方面展開研究[5]。陽艾利等[6]研究表明，Gumbel Copula函數可用于降雨量和徑流量聯合分布的擬合；陳心池等[7]研究表明，Frank Copula函數可用于漢江中上游流域極端降雨洪水三維聯合分布擬合。彭楊等[8]基于多維Copula函數進行日含沙量過程隨機模擬，模擬結果表明，該法能保持實測日含沙量的主要統計特性，模擬序列與實測序列的統計特征值吻合良好。李帆等[9]基于Archimedean Copula函數對淮河月徑流進行隨機模擬，研究結果表明，Archimedean Copula函數結合Gibbs采樣可有效建立相鄰月份間的相關性結構。

相對于傳統的單一站點預報，通過協整理論建立多站點時空關系，并運用非線性關系進行聯合預報的多站點同步預報方式，具有更為突出的預報效果。吳金冉等[10]以青海省海東和果洛月降水量數據為例，得出結論：多站點同步預報模型，相對于傳統單站點預報模型具有更為突出的預報效果。

本文以江蘇省宿遷市古黃河流域典型站點實測降雨量為研究對象，基于Copula函數[11]進行3個站點兩兩間月降雨量聯合分布擬合，并結合Gibbs抽樣法[12]，實現月降雨量隨機模擬。通過比較實測和模擬月降雨量，對模擬結果的有效性進行驗證，為多站點間的中長期概率水文預報提供參考。

1 數據與方法

1.1 數據來源

宿遷市位于江蘇省北部。古黃河是宿遷市主要河流之一，其橫穿宿遷市城區中部，全長114.3 km。流域面積約296.9 km2。古黃河流域屬于暖溫帶季風性氣候，四季分明，年際降雨量變化較大，且年內分布不均，易形成春旱、夏澇、秋冬干天氣。多年平均降水量為913.1 mm。多年汛期平均降雨量611.5 mm，約占多年年平均降雨量的67.0%，最大汛期降雨量1 256.0 mm(2003年)，最小降雨量283.1 mm(1966年)。本文選取分布于古黃河流域上、中、下游的宿遷閘站、劉老澗站、新袁閘站3個代表站(圖1)39年(1980—2018年)月實測降雨量資料作為研究對象。

圖1 研究區域及站點分布

1.2 二維Copula函數理論

Copula函數的種類較多，常用的有橢圓Copula函數中的Gaussian Copula函數和Archimedean Copula函數中的Gumbel Copula、Clayton Copula、Frank Copula函數[13]。其中Archimedean Copula函數因其適用性強，構造簡單、計算簡便，應用十分廣泛。

常用二維Archimedean Copula函數如下：

Gumbel CopulaC(u,v)=e-[(-lnu)θ+(-lnv)θ]1/θ

(1)

Clayton CopulaC(u,v)=(u-θ+v-θ-1)-1/θ

(2)

(3)

Copula 函數的參數估計方法有很多[14]，目前較常用的有非參數法和兩階段估計法[15]，本文采用非參數法進行參數求解。

非參數法通過構建Kendall秩相關系數τ與Copula函數參數θ間相關關系來求解參數：

(4)

(5)

(6)

1.3 分布函數優選

本文選取P-III型分布、指數分布和Gamma分布作為待選邊緣分布。通過計算各邊緣分布均方根誤差來選取最優分布，并通過 Kolmogorov-Simirnov檢驗(K-S檢驗)對樣本進行擬合優度檢驗。再選取Gumbel Copula、Clayton Copula、Frank Copula函數為備選函數建立模型，通過Kendall秩相關系數τ與式(4)—(6)，計算相關參數。采用赤池信息法(Akaike Information Criterion，AIC)[16]優選Copula 函數模型，AIC最小為最優。

2 流域內站點間月降雨量聯合分布模型建立

2.1 各站降雨量邊緣分布確定

目前，中國水文計算分析中常常假定降雨、徑流服從的最優分布為P-III分布。但實際情況中，P-III型分布并非具有普適性。在Copula函數中，選取擬合優度最好的邊緣分布，是決定聯合分布模型擬合精度的關鍵因素之一。因此，選取P-III型分布、指數分布和Gamma分布作為待選邊緣分布，通過計算各邊緣分布均方根誤差來選取最優分布(表1)。

表1 邊緣分布擬合優度檢驗結果

由表1可知，根據RMSE最小準則，3個站點1—12月降雨變化P-III分布擬合的RMSE值較Gamma分布和Expoerntial分布小，因此，選取P-III分布作為宿遷閘站、劉老澗站、新袁閘站月尺度降雨變化的邊緣分布。其中，K-S檢驗值均大于0.05，統計值均小于臨界值，故認為通過檢驗。

2.2 Copula函數的選擇

本文以Archimedean Copula函數中的Gumbel Copula、Clayton Copula、Frank Copula作為備選函數建立聯合概率分布函數。由Gumbel Copula邏輯模型的結構可知，其適用于兩變量之間存在正相關關系的序列；Clayton Copula函數一般適用于有正相關關系的隨機變量；Frank Copula函數適用于正負相關的隨機變量。根據3個站點降雨變化相關性的大小來初步選取Copula函數。為此，首先分析各個相鄰站點間降雨量相關性大小。

由表2可知，上、中、下游3個站點的兩兩相關性較好，降雨相關系數最大達0.99，說明3個站降雨變化總體上存在較強正相關。

表2 相關系數計算結果

2.3 Copula函數參數估計及擬合優度檢驗

二維Copula函數的參數根據式(4)—(6)進行估計。為確定最優Copula函數，需對初選的3個Copula函數進行擬合優度檢驗。對3個站點不同時段降雨變化的觀測值進行兩兩組合，(xi,yi)按x升序排列，得到新數據組(x1,y1)，(x2,y2)，…，(xn,yn)，二維經驗聯合概率計算公式為：

Fexp(xi,yi)=p(X≤xi,Y≤yi)=

(7)

式中Fexp——經驗聯合概率分布值；ng,k——所有滿足X≤xi、Y≤yi的聯合觀測值的數目；n——序列長度。

其中5、9月降雨變化經驗聯合頻率與理論聯合頻率的一致性檢驗見圖2—5。從圖可知，Clayton Copula、Frank Copula、Gumbel Copula 函數構造的理論頻率和分布模型對經驗聯合頻率的一致性較強，線性相關性在0.98以上。需要通過其他方法對Copula函數進行優選。

為選取最優Copula函數，采用AIC最小準則進行擬優選(表3)。通過分析，Gumbel Copula函數普遍優于Clayton Copula函數、Frank Copula函數。

a)Clayton Copula

a)Clayton Copula

a)Clayton Copula

a)Clayton Copula

表3 AIC準則擬合優度檢驗結果

3 Copula函數隨機模擬

隨機模擬常采用的模擬方法主要有直接抽樣法(一般適用于一維)、拒絕抽樣法、重要性抽樣法、MCMC抽樣方法(分為Metropolis-Hasting算法與Gibbs采樣算法)。

無論是拒絕抽樣還是重要性采樣，都是屬于獨立采樣，即樣本與樣本之間是獨立無關的，這樣的采樣效率比較低；MCMC方法是關聯采樣，即下一個樣本與這個樣本有關系，從而使得采樣效率高。

Metropolis-Hasting算法的使用需要滿足細致平衡的要求，而且這個算法有一個缺點，就是抽樣的效率不高，有些樣本會被舍棄掉，從而產生了Gibbs算法。Gibbs算法是用條件分布的抽樣來替代全概率分布的抽樣，而且舍棄的樣本較少，效率較高，且能夠較好地保留原有數據的屬性特征。

Copula函數模擬月降雨量，采用條件概率的反函數，結合Gibbs抽樣法對3個水文站的月降雨量進行隨機模擬。根據選定的Copula函數，構造2個站之間的聯合分布，利用Copula函數計算條件分布。

由表3結果可得，本次所選最優Copula函數中無Clayton Copula，故下文僅對Frank Copula、Gumbel Copula進行分析。根據式(1)—(3)可以得到Franke Copula和Gumbel Copula函數的條件分布函數，分別為：

F(v|u)=?C(v,u)/?u=e-uθ(e-vθ-1)/

[(e-vθ-1)(e-uθ-1)+e-θ-1]

(8)

F(v|u)=?C(u,v)/?u=

(9)

利用上述條件分布，結合Gibbs采樣方法對月降雨量進行隨機模擬。Gibbs采樣方法的步驟如下。

步驟一生成隨機數a0、a1∈(0,1);

步驟二令P(X1,1,1≤x1,1|X2,1,1=a0)=?C(x1,1,1,x2,1，1)/?x2,1,1=a1,其中Xi,j,k表示第k年j月i站的月降雨量，將其代入式(1)，求出v即為x1,1,1；

步驟三生成2個隨機數a2、a3∈(0,1),求解方程組：

(10)

即得第1年第1月3個站點的月降雨量x1,1,1、x2,1,1、x3,1,1。

重復步驟三500次，可得500年的模擬月降雨量(圖6、7，黑點表示實測值，灰點表示模擬值)。

a)宿遷閘站-劉老澗站

a)宿遷閘站-劉老澗站

4 結果分析

由圖7可得，3站點間模擬值與實測值的走向基本一致，模擬點與實測點也較為集中，說明擬合值和實測值差別較小，擬合效果較好。

為進一步分析模擬值與實測值間的擬合效果，進行模擬值與實測值的均值、均方差、變差系數的誤差分析，結果見表4。由表4可知，1月宿遷閘站與劉老澗站、3月劉老澗站、4月宿遷閘站的均方差、變差系數相對誤差較高，經過分析可知，原數據存在降雨量極大值，導致邊緣分布擬合誤差較大。模擬值與實測值的均值、均方差、變差系數整體相對誤差較小，基本處于10%以下，模擬效果較好。

表4 研究區模擬值與實測值誤差分析

相鄰站點間相關系數的誤差分析結果見表5。

續表4 研究區模擬值與實測值誤差分析

表5 研究區相鄰站點間相關系數模擬值與實測值的誤差

續表5 研究區相鄰站點間相關系數模擬值與實測值的誤差

由表5可知，各站點間Kendall相關系數、Spearman相關系數相對誤差均較小，最大為8.93%。說明本次模擬基本保留了原有觀測值的統計特征。

5 結論

本文通過Copula函數法對上游宿遷閘站、中游劉老澗站、下游新袁閘站3站的月降雨量進行了隨機模擬，結果表明，該3站月降雨量存在較強的正相關性，二維Copula函數結合Gibbs采樣法可有效建立同流域相鄰站點間的相關性結構，模擬結果能夠較好地保留原有觀測值的統計特征。該法簡單易懂，可為其他水文要素的隨機模擬與中長期概率水文預報提供參考。