北方干旱半干旱區洪水極值事件多變量頻率分析

曹 丹

(遼寧省葫蘆島水文局，遼寧 葫蘆島 125000)

在水文分析計算中，常常根據某一類極值事件，例如洪水、暴雨，單變量，對該極值事件中包含的某一種特征屬性，如洪峰、洪量、洪水持續時間進行頻率分析以得到該類極值事件下某一頻率對應的設計值[1]。然而在實際中，一種屬性往往無法準確地描述該類極值事件的具體特征，也無法構建該類極值事件的某一特征和其他事件的某一屬性值之間的相互關系。一些特征變量之間的聯合分布無法采用單變量頻率分析方法進行全面描述，變量之間的相關結果也無法刻畫。因此，要解決此項問題需要通過多變量水文頻率分析的方式進行[2]。Copula函數是近年來常用于多變量頻率分析的一種方法[3- 12]，與其他多變量頻率分析方法相比，Copula函數對于邊緣獨立分布的隨機變量進行相關性結構的反映，分為2個獨立的部分來對變量聯合分布進行計算，從而構建多個相關結構的變量邊緣分布概率函數，通過函數變化形成多個變量之間的聯合分布，由于在邊緣分布函數里對變量所有特征進行包含，因此在函數變換過程中信息不會有所損失。通過Copula函數可以構建洪峰流量、洪水歷時、洪量之間的聯合分布來全面地描述洪水，也可以將洪水事件與降水事件結合，通過考慮氣象要素與洪水特征之間的相互關系，進而分析氣象要素對洪水極值事件的影響[13]。從而對某一區域的洪水頻率分析提供更全面、更科學的設計值[14]。遼西地區近些年來受氣候變化影響，洪水極值發生不同程度的變化，亟需對其洪水極值和降雨條件聯合分布概率進行探討，從而為應對氣候變化影響的區域洪水影響制定相對應的科學措施。為此本文選取東白城子以上流域為研究實例，結合區域實測洪水和降水數據，基于二維Copula函數建立洪水極值和降水的聯合分布概率，從而確定區域洪水極值在不同降水條件下的合概率分布。

1 研究方法

Copula函數的種類較多，其中Gumbel-Hougaard Copula函數的上尾和下尾相關系數適合對水文變量進行描述，其函數主要原理由于參考文獻較多，本文就不做過多介紹，分別采用Pearson相關系數等5種方法對其變量間的相依性度進行計算，Copula函數概率分析時候需要設{(x1，y1)，(x2，y2)，……(xn，yn)}為分析變量x、y的連續觀測樣本，則各相關系數的計算方程分別為：

(1)

(2)

(3)

(4)

除上述3種相關系數外，Chi圖和K圖也是常見的用于描述變量間相依性的方法，Chi圖的計算原理如下：

(5)

(6)

(7)

則有：

(8)

(9)

(10)

(11)

式中，

(12)

2個變量X、Y在K圖中具有一定的相依度，Win和K0(ω)關系曲線均會對(Win，H(i))的所有數據在帶狀區域內進行圈圍，2個變量之間的越高的相依度則其關系曲線越為接近。

2 研究區域概況

東白城子水文站以上流域位于繞陽河流域內，繞陽河源出阜新蒙古族自治縣扎蘭營子鄉察哈爾山，于盤山縣東郭鎮萬金灘注入遼河，河流總長度326km，流域面積10438.20km2，其中干流河段河道長度48.19km。流域內山丘漫崗區占總面積的43.30%，沖積平原區占56.70%，平均坡降1.72%。流域降水量和蒸發量多年均值分別為485、1750mm，降雨受季風氣候影響時空分布差異較大，從南到北降水逐步遞減，降水主要集集中在汛期，占全年降水量的比重超過80%。低山丘陵是繞陽河流域的主要地形，土壤產水能力較低，暴雨洪水主要集中在汛期。繞陽河干流自河源以下共有4個水文站，分別是韓家杖子站、東白城子站、四家子站和杜家站。東白城子水文站以上流域為繞陽河上游，集水面積2070km2。東白城子流域內共有10個配套雨量站。

3 分析結果

3.1 變量序列樣本選取

選取的水文站點為東白城子流域出流站東白城子站。基于東白城子水文函及流域內10個雨量站時間序列的交集，選取的時間序列范圍為1967—2018年。采用年最大值法(AM)選取東白城子站年最大洪峰序列作為洪水極值事件的洪峰序列，在降水樣本選取時，充分考慮降水和洪水之間的相關性以及兩者在物理成因上的相互而選取最具代表性的降水序列。基于此，本節選用前期累計面降水量序列作為降水序列。采用泰森多邊形法，基于東白城子流域10個雨量站點的點降水量資料，計算東白城子流域的面降水量序列。隨后，從最大洪峰發生的時間開始，從當天往前推算Nd的累積降水量，通過計算不同前推時間N下，累積降水量序列和年最大洪峰序列之間的相關系數，選取相關系數對應的Nd累積降水量序列進行下一部分的工作。前期累積降水量(N=1，2，……)和年最大洪峰流量序列之間的相關系數如圖1所示。

圖1 前期累積面雨量序列與年最大洪峰序列的相關性

可以清晰看出，當N=18d時，前期累積降水量與年最大洪峰序列之間的相關性最高，其對應的相關系數為0.792。N=18d對應的相關性如圖1所示。為此最大洪峰年序列和降水量前期累積量對應的N值為18d，2個變量之間的相應性結構采用Copula函數進行描述。

3.2 多變量邊緣分布的確定

采用極大似然法估計確定東白城子站1960—2018年年最大洪峰流量序列、年最大洪量序列、年最長洪水歷時序列10種頻率分布下的參數，并采用擬合度比較方法確定各個序列的最優分布。東白城子站、韓家杖子站各分布的經驗頻率與理論頻率的擬合度及最優分布選擇結果見表1。

表1 東白城子站洪峰、洪量、歷時最優分布選擇

從計算結果中可以看出，在東白城子站，1960—2018年年最大洪峰序列、年最大洪量序列、年最長洪水歷時序列的最優分布分別是P-Ⅲ型分布、GPD分布。基于此，選擇各變量對應的最優分布作為構建Copula函數的邊緣分布。

3.3 序列的相依性度量計算

年最大洪峰和洪量序列的分析年份為1960—2018年，分別對洪峰和洪量、洪峰與歷時以及洪量與歷時的Pearson相關系數γ，Spearman相關系數ρ和Kendall相關系數τ進行兩兩序列計算，結果見表2。對于相依性較好的洪峰與洪量序列，可以采用Copula函數描述兩者之間的相依性結構，對于相關關系不好的序列，不建議采用Copula函數進行分析。為此還采用Chi圖、K圖進一步驗證東白城子站年最大洪峰流量、年最大洪量之間的相關性，如圖2所示。

表2 東白城子站兩變量間的相依性度量

圖2 東白城子站年最大洪峰洪量序列Chi和k圖

從計算結果中可以看出，對于東白城子站，洪峰、洪量之間的相依性較好，不論是Pearson相關系數γ，Spearman相關系數ρ和Kendall相關系數τ均達到了0.6以上，而洪峰與洪水歷時，洪量與洪水歷時之間的相依性并不好，不論是Pearson相關系數γ，Spearman相關系數ρ和Kendall相關系數τ，均沒有超過0.2。從圖2中也可以看出，東白城子站年最大洪峰序列與年最大洪量序列之間具有良好的相依性，在其Chi圖中，除個別點外，其余所有點據均落在虛線置信區間外。而在K圖中，除個別點據外，Win和K0(ω)的關系曲線和對角線所圍區域內均包含剩余點據，關系曲線和點據具有較好的擬合度，由此說明東白城子站年最大洪峰序列與年最大洪量序列之間具有較強的相依性。此外，從Chi圖中還可以發現，除了一個點外，所有的點的χi均大于0，由此也進一步說明了東白城子站的年最大洪峰序列與年最大洪量序列之間存在正象限相關關系。

3.4 擬合度檢驗結果

分別采用Gumbel Copula，Clayton Copula，Frank Copula 3種Copula函數，建立東白城子站年最大洪峰、洪量序列之間的聯合分布函數，分別采用2種方法對不同Copula函數進行參數估計，最后進行擬合優度檢驗并選擇擬合效果最優的Copula函數。根據RMSE的計算公式，RMSE值越小代表擬合優度越高，以此為標準選擇最優的Copula函數及最優的參數估計方法。東白城子站年最大洪峰洪量序列Copula參數計算結果見表3。

表3 東白城子站年最大洪峰洪量序列的Copula函數參數估計及擬合優度檢驗

從上表中可以看出，采用矩估計法(MOM法)對東白城子站年最大洪峰序列與年最大洪量序列之間不同Copula函數進行參數估計時，擬合效果最佳的為Gumbel Copula函數，擬合效果相對較差的是Frank Copula函數和Clayton Copula函數。在該站采用偽極大似然法(MPL法)估計參數時，采用Clayton Copula函數的擬合效果最好，Frank Copula函數擬合效果最低，Gumbel Copula函數擬合效果好于Frank Copula函數，因此在東白城子站采用MPL方法得到的Copula函數的擬合效果均更佳。此外，在東白城子站，通過MPL方法估計Clayton Copula函數參數得到的擬合效果最佳，由此說明Clayton Copula函數為東白城子站年最大洪峰、洪量序列的最優擬合分布。因此選用Clayton Copula函數，用于構建東白城子站年最大洪峰序列、年最大洪量序列之間的聯合分布函數，且參數估計方法均選用MPL方法。

3.5 多變量模型分析

基于東白城子站、韓家杖子站的年最大洪峰、洪量之間的最優Copula函數，繪制了相應的聯合重現期等值線圖(重現期分別為5、10、20、50、100年)，其結果如圖3所示。

圖3 東白城子站洪峰洪量兩變量聯合重現期

年最大洪量/洪峰流量無論何種重現期，隨著年最大洪峰流量/洪量的增加其發生的條件概率呈現遞增變化，即說明東白城子站年最大洪峰與年最大洪量之間呈現明顯的正相關關系，在洪峰流量較大的情況下，發生較大洪量洪水的可能性也較大，這也與實際事實相符合。此外，從上圖中也可以看出，在同一洪峰流量/洪量下，重現期較小的情況下對應的洪量/洪峰流量發生的累積概率要明顯大于重現期較大情況下對應的洪量/洪峰流量發生的累積概率，這一點與實際情況也是相符的，因為在相同的洪峰流量/洪量條件下，重現期小的洪水比重現期大的洪水更易發生。

4 結論

(1)采用矩估計法(MOM法)估計北方干旱半干旱區域年最大洪峰序列與年最大洪量序列之間不同Copula函數的參數時，擬合效果最好的是Gumbel Copula函數，擬合效果相對較差的是Clayton Copula函數和Frank Copula函數。

(2)對于北方干旱半干旱區而言，洪峰、洪量之間的相依性較好，Pearson相關系數γ、Spearman相關系數ρ和Kendall相關系數τ均達到了0.6以上，而洪峰與洪水歷時，洪量與洪水歷時之間的相依性較低，Pearson相關系數γ、Spearman相關系數ρ和Kendall相關系數τ均低于0.2。

(3)本文在變量樣本選取時只考慮了降水和洪水極值之間的相關性，在后續的研究中還應增加其他氣象要素進行洪水極值樣本的綜合選取，增加樣本序列的覆蓋程度。