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

曹 丹

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

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

1 研究方法

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

(1)

(2)

(3)

(4)

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

(5)

(6)

(7)

則有：

(8)

(9)

(10)

(11)

式中，

(12)

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

2 研究區(qū)域概況

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

3 分析結(jié)果

3.1 變量序列樣本選取

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

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

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

3.2 多變量邊緣分布的確定

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

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

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

3.3 序列的相依性度量計(jì)算

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

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

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

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

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

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

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

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

3.5 多變量模型分析

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

圖3 東白城子站洪峰洪量?jī)勺兞柯?lián)合重現(xiàn)期

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

4 結(jié)論

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

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

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