曹 丹

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

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

1 研究方法

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

(1)

(2)

(3)

(4)

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

(5)

(6)

(7)

則有：

(8)

(9)

(10)

(11)

式中，

(12)

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

2 研究區(qū)域概況

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

3 分析結(jié)果

3.1 變量序列樣本選取

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

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

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

3.2 多變量邊緣分布的確定

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

表1 東白城子站洪峰、洪量、歷時最優(yōu)分布選擇

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

3.3 序列的相依性度量計算

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

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

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

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

3.4 擬合度檢驗結(jié)果

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

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

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

3.5 多變量模型分析

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

圖3 東白城子站洪峰洪量兩變量聯(lián)合重現(xiàn)期

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

4 結(jié)論

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

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

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