基于Copula函數(shù)的設(shè)計(jì)暴雨雨型研究

葉姍姍， 葉興成， 王以超， 朱程亮， 劉 俊

(1.河海大學(xué) 水文水資源學(xué)院， 江蘇 南京 210098； 2.宿遷市水務(wù)局， 江蘇 宿遷 223800)

1 研究背景

隨著全球氣候的不斷變化，城市化建設(shè)進(jìn)程加快和原有排水防澇系統(tǒng)設(shè)計(jì)標(biāo)準(zhǔn)偏低等因素，宿遷市常出現(xiàn)較嚴(yán)重的洪澇災(zāi)害，因此科學(xué)確定城市排澇設(shè)施的建設(shè)標(biāo)準(zhǔn)已成為解決城市暴雨內(nèi)澇的重要任務(wù)之一。設(shè)計(jì)雨型是制定排水系統(tǒng)設(shè)計(jì)標(biāo)準(zhǔn)的重要依據(jù)和前提，能為城市排水系統(tǒng)的規(guī)劃和管理、暴雨洪水分析和流域調(diào)洪計(jì)算提供科學(xué)依據(jù)[1]。國內(nèi)的雨型研究開展較晚，以往常用同倍比和同頻率雨型，后又引入Keifer和Chu雨型、Pilgrim & Cordery雨型、huff雨型和Yen & Chow雨型[2]，岑國平等[3]對(duì)此進(jìn)行了大量的研究。國外的雨型設(shè)計(jì)除上述的方法外還有概率密度函數(shù)法、隨機(jī)模型法和離散法等。隨機(jī)模型法分為單站點(diǎn)模型(如Newman-scott模型和Bartlett-Lewis模型)、多站點(diǎn)模型(如馬爾可夫鏈模型)和時(shí)空模型。離散法有級(jí)聯(lián)過程法、混沌過程法和人工神經(jīng)網(wǎng)絡(luò)法等[4]。

近年來國內(nèi)的雨型研究進(jìn)展緩慢，還是以芝加哥雨型和PC雨型為研究重點(diǎn)，和國外的雨型研究成果相差甚遠(yuǎn)，進(jìn)一步深入研究和推進(jìn)設(shè)計(jì)雨型工作具有重大意義。因此本文以宿遷市1981-2016年36 a的實(shí)測逐時(shí)降雨資料為例，在原有設(shè)計(jì)雨型方法的基礎(chǔ)上進(jìn)行創(chuàng)新，應(yīng)用同頻率法SCS法、暴雨衰減指數(shù)法和PDF法來建立宿遷設(shè)計(jì)雨型，并用Copula函數(shù)的風(fēng)險(xiǎn)聯(lián)合概率將其與同頻率雨型進(jìn)行對(duì)比，分析這3種方法應(yīng)用于宿遷市設(shè)計(jì)暴雨雨型研究的合理性。

2 設(shè)計(jì)雨型的研究方法

宿遷地處亞熱帶和暖溫帶的過渡地區(qū)，主要受北方西風(fēng)槽、冷渦和熱帶臺(tái)風(fēng)，及當(dāng)?shù)厣傻慕辞凶兙€和氣旋等天氣系統(tǒng)影響，暴雨主要發(fā)生在7-9月份，具有時(shí)長、量大、面廣的特點(diǎn)，其中特大暴雨會(huì)對(duì)其城市防洪排澇系統(tǒng)帶來巨大挑戰(zhàn)。宿遷閘站位于宿遷中心城區(qū)邊緣，離中心城區(qū)最近，是水文部門的國家基礎(chǔ)站，其可靠性、一致性和代表性良好。

選取宿遷閘站觀測數(shù)據(jù)，采用年最大值法選取1981-2016年的實(shí)測24 h降雨資料，運(yùn)用下面的4種設(shè)計(jì)雨型方法推求宿遷市設(shè)計(jì)暴雨雨型。

2.1 同頻率法雨型

同頻率法即選取當(dāng)?shù)貙?shí)測的典型暴雨，對(duì)其降雨過程進(jìn)行同頻率分時(shí)段縮放[5]。放大倍比系數(shù)公式：

(1)

式中：X設(shè)計(jì)為不同重現(xiàn)期的設(shè)計(jì)降雨量，mm；X典型為典型暴雨降雨量，mm。

宿遷市2010年9月7日遭遇特大暴雨，日雨量為歷史第二位，達(dá)35年一遇，引起巨大的洪澇災(zāi)害。該次暴雨雨型為主副型，雨量集中，雨峰偏后，在水文氣象特征上具有一定代表性，故選取該特大暴雨過程作為典型暴雨。

采用皮爾遜III型目估適線法[6]計(jì)算宿遷市各重現(xiàn)期各時(shí)段設(shè)計(jì)降雨量如表1。

表1 不同重現(xiàn)期下各歷時(shí)設(shè)計(jì)降雨量 mm

典型暴雨的最大1、3、6、12和24 h的時(shí)段降雨量分別為65.9、147.2、221.1、251.1和256.7 mm，根據(jù)公式(1)可計(jì)算出各重現(xiàn)期的各時(shí)段放大系數(shù)，通過縮放典型暴雨過程即得24 h的設(shè)計(jì)雨型如圖1所示。

圖1 宿遷市24 h設(shè)計(jì)雨型(同頻率法)

2.2 SCS雨型

SCS雨型是由美國水土局于20世紀(jì)80年代建立的，具體分為4種類型。其中SCS-Ⅱ代表由夏季陣雨造成的降雨歷時(shí)短但強(qiáng)度高的降雨地區(qū)，與宿遷市的降雨特征較相似[7]。SCS-Ⅱ雨型的各時(shí)段百分比累積數(shù)據(jù)見表2。

表2 SCS-Ⅱ雨型時(shí)程分布數(shù)據(jù)

由表2各時(shí)段累積比例計(jì)算降雨時(shí)程分配比例，根據(jù)宿遷市不同重現(xiàn)期的24h設(shè)計(jì)雨量可推求出其SCS-Ⅱ設(shè)計(jì)雨型如圖2所示。

圖2 宿遷市24 h設(shè)計(jì)雨型(SCS-Ⅱ法)

2.3 暴雨衰減指數(shù)法雨型

宿遷閘站有詳細(xì)的1、3、6和24 h的降雨系列且屬于中小流域，可引入暴雨衰減指數(shù)來推求設(shè)計(jì)雨型[8]。暴雨衰減指數(shù)越大表示降雨越集中。由小流域暴雨公式可知，第k時(shí)段暴雨量hk為：

(2)

式中：hk為匯流時(shí)間(t)內(nèi)的暴雨量，mm；ik為匯流時(shí)間內(nèi)的暴雨強(qiáng)度，mm/h；H24p為24 h的暴雨量，mm；Bk為暴雨分配系數(shù)；np為暴雨衰減系數(shù)。

當(dāng)Bk=1時(shí)，反推算匯流時(shí)間區(qū)間為(t1，t2)的暴雨衰減系數(shù)為：

(3)

式中：t1，t2為匯流時(shí)間，h；ht1，ht2為對(duì)應(yīng)匯流時(shí)間的設(shè)計(jì)降雨量，mm。

由表1宿遷市的設(shè)計(jì)雨量，根據(jù)式(3)可以計(jì)算其暴雨衰減指數(shù)如表3。

表3 不同時(shí)間區(qū)間的np值

已知np的前提下，各時(shí)段設(shè)計(jì)雨量可根據(jù)以下公式計(jì)算[9]：

(1)ti=1～3 h

Hi=H1·ti1-np(1,3)或

Hi=H3·(ti/3)1-np(1,3)

(4)

(2)ti=3～6 h

Hi=H3·(ti/3)1-np(3,6)或

Hi=H6·(ti/6)1-np(3,6)

(5)

(3)ti=6～24 h

Hi=H6·(ti/6)1-np(6,24)或

Hi=H24·(ti/24)1-np(6,24)

(6)

根據(jù)《浙江省可能最大暴雨圖集》可知，將計(jì)算得出的各時(shí)段設(shè)計(jì)雨量的最大值放在峰值位置，次大降雨量放在最大值左邊，其余時(shí)段雨量從大到小排序，奇數(shù)放在左邊，偶數(shù)放在右邊，當(dāng)一邊排滿時(shí)其余的按從大到小全放到另一邊。由宿遷市36 a的24 h的降雨過程計(jì)算得雨峰時(shí)刻平均值為第10 h，由此得到設(shè)計(jì)雨型如圖3所示。

圖3 宿遷市24 h設(shè)計(jì)雨型(暴雨衰減指數(shù)法)

2.4 PDF法雨型

2.4.1 PDF法雨型簡介及參數(shù)推求 研究得知累積降雨曲線與累積分布函數(shù)(CDF)形狀相似，而累積分布函數(shù)(CDF)的一次微分為概率密度函數(shù)(PDF)，因此可用PDF來設(shè)計(jì)雨型。而且國外已經(jīng)初步研究了用PDF來構(gòu)建雨型，結(jié)果發(fā)現(xiàn)Gumbel、Lognormal和Weibull等概率密度函數(shù)可很好的反映降雨過程特征，適用于推求設(shè)計(jì)雨型[4]。其具體型式如下：

(1)Gumbel分布

(7)

式中：α與β為參數(shù)。

(2)Lognormal分布

(8)

式中：η與σ為參數(shù)。

(3)Weibull分布

(9)

式中：γ與α為參數(shù)。

以Gumbel分布為例，公式如下：

Tp=α

(10)

(11)

(12)

參數(shù)α，β要同時(shí)滿足上述的3個(gè)目標(biāo)函數(shù)，顯然這是一個(gè)多元非線性方程組的求解問題，不能直接解出數(shù)值，需利用最優(yōu)化算法進(jìn)行求解。因此本文采用絕對(duì)值求和法，將多目標(biāo)函數(shù)問題轉(zhuǎn)化成單目標(biāo)函數(shù)最優(yōu)化問題，只要求解出絕對(duì)值之和最小時(shí)的參數(shù)即可，參數(shù)結(jié)果應(yīng)滿足下式：

(13)

2.4.2 粒子群優(yōu)化算法 PDF參數(shù)求解問題轉(zhuǎn)化成了單目標(biāo)函數(shù)minZ的最優(yōu)化求解問題。粒子群算法(PSO)算法簡單，求解速度快，適合求解實(shí)數(shù)問題，因此本文采用PSO對(duì)函數(shù)minZ進(jìn)行最優(yōu)化求解，解出Z最小時(shí)對(duì)應(yīng)的參數(shù)結(jié)果。

粒子群優(yōu)化算法采用速度-位置模型，每個(gè)粒子代表著參數(shù)解空間的一個(gè)，解的好壞程度由適應(yīng)函數(shù)的結(jié)果來決定[10]。本文的粒子群算法的學(xué)習(xí)因子取c1=c2=2。為了使粒子不遠(yuǎn)離解空間，粒子移動(dòng)速度應(yīng)控制在[-vmax，vmax]，一般取vmax=k·vmax，本文取k=0.5。粒子在解空間不斷更新迭代直到達(dá)到規(guī)定的迭代次數(shù)或達(dá)到規(guī)定的誤差標(biāo)準(zhǔn)為止[11]。

應(yīng)用matlab進(jìn)行粒子群算法編程，以宿遷市2000年的數(shù)據(jù)為例，計(jì)算出Tp、Qp、Qd的輸入值，分別用Gumbel、Lognormal和Weibull概率密度函數(shù)來擬合2000年實(shí)測降雨過程，求解minZ最小，由求解的參數(shù)得擬合圖如圖4所示。

采用粒子群最優(yōu)化算法分別擬合出每年的Gumbel、Lognormal和Weibull的分布曲線。以宿遷市1981-2016的實(shí)測降雨數(shù)據(jù)為基礎(chǔ)進(jìn)行函數(shù)間的擬合優(yōu)度檢驗(yàn)，選擇出與實(shí)測數(shù)據(jù)擬合最好的函數(shù)。本文擬合優(yōu)度檢驗(yàn)采用3種方法[12]：

(1)擬優(yōu)平方和準(zhǔn)則法

(14)

(2)擬優(yōu)絕對(duì)值準(zhǔn)則法

(15)

(3)概率點(diǎn)據(jù)相關(guān)系數(shù)檢驗(yàn)法

(16)

式中：x(i)為函數(shù)擬合值；y(i)為實(shí)測樣本值；xm為擬合值的均值；ym為樣本值的均值。

由上述公式計(jì)算出1981-2016年優(yōu)度檢驗(yàn)均值，結(jié)果見表4。擬合優(yōu)度檢驗(yàn)中RMSE和MAE越小越優(yōu)，PPCC越大越優(yōu)。對(duì)于宿遷市實(shí)測數(shù)據(jù)來說，Gumbel分布是最優(yōu)分布，可作為其設(shè)計(jì)雨型。

表4 擬合優(yōu)度檢驗(yàn)結(jié)果

2.4.3 PDF法設(shè)計(jì)雨型結(jié)果 選用Gumbel分布推求宿遷市設(shè)計(jì)雨型，由宿遷市36 a的24 h的降雨過程得雨峰時(shí)刻平均值為第10 h，即輸入值Tp=10。由宿遷市不同重現(xiàn)期下24 h設(shè)計(jì)雨量和最大1 h設(shè)計(jì)雨量(如表2)得輸入值Qp、Qd。運(yùn)用粒子群算法計(jì)算Gumbel分布參數(shù)，得設(shè)計(jì)雨型如圖5所示。

3 基于Copula函數(shù)設(shè)計(jì)雨型的比較分析

不同設(shè)計(jì)雨型的優(yōu)劣沒有固定的方法來判別，本文嘗試?yán)肅opula函數(shù)，從設(shè)計(jì)降雨總量分別與最大1h雨量、最大3h雨量、最大6h雨量間的風(fēng)險(xiǎn)聯(lián)合分布概率入手，對(duì)上述4種雨型進(jìn)行對(duì)比分析。

Copula函數(shù)是定義在[0，1]區(qū)間均勻分布的聯(lián)合概率分布函數(shù)[13]。由宿遷市36 a實(shí)測數(shù)據(jù)得出變量間為正相關(guān)，故本文選取了Gumbel-Hougaard Copula函數(shù)來分析設(shè)計(jì)雨型。GH Copula函數(shù)的表達(dá)式為：

(17)

θ∈[1,∞]

式中：θ為參數(shù)。

采用皮爾遜III型來計(jì)算邊緣分布函數(shù)，由矩法得出P-III參數(shù)值見表5，通過目估適線法可推求出不同降雨量對(duì)應(yīng)的頻率值(即u、v的值)。

(18)

式中：(xi，yi)、(xj，yj)分別為兩組獨(dú)立同分布的觀測數(shù)據(jù)；sign為符號(hào)函數(shù)。

表5 P-III邊緣分布函數(shù)的參數(shù)值

計(jì)算得24 h降雨量與最大1、3、6 h降雨量的Kendall秩相關(guān)系數(shù)τ分別為0.5016、0.5952和0.6302，則對(duì)應(yīng)GH Copula的參數(shù)θ分別為2.01、2.47和2.7。

以24 h降雨量與最大1 h降雨量的Copula概率分布關(guān)系為例來詳細(xì)介紹。采用Genest-Rivest法驗(yàn)證理論聯(lián)合分布函數(shù)與經(jīng)驗(yàn)聯(lián)合分布函數(shù)的擬合程度，結(jié)果如圖6所示。各點(diǎn)均勻的分布在45°線左右，則GH Copula具有良好的擬合度[15-16]。

本文定義風(fēng)險(xiǎn)率為滿足24 h設(shè)計(jì)降雨量條件下，最大1h(或3、 6 h)降雨量超出對(duì)應(yīng)設(shè)計(jì)時(shí)段降雨量發(fā)生的條件概率[5]，計(jì)算公式如下式：

(1)

式中：R為24 h設(shè)計(jì)降雨量，mm；Ri為最大1h(或3、 6 h)時(shí)段設(shè)計(jì)降雨量，mm；F(R)為設(shè)計(jì)降雨量r的分布函數(shù)；F(R，Ri)為(R，Ri)的聯(lián)合分布函數(shù)。

圖4 2000年3個(gè)分布函數(shù)的擬合曲線圖

根據(jù)GH Copula函數(shù)公式可計(jì)算出宿遷閘站24 h降雨量與最大1 h降雨量的不同組合下的聯(lián)合分布概率F(R，Ri)，以此可計(jì)算出其風(fēng)險(xiǎn)聯(lián)合分布概率[17]，結(jié)果如圖7～8所示。

圖5 宿遷市24 h設(shè)計(jì)雨型(PDF法雨型)

圖6 經(jīng)驗(yàn)聯(lián)合分布函數(shù)與理論聯(lián)合分布函數(shù)

圖7 24 h降雨量與最大1h降雨量的聯(lián)合頻率分布圖

圖8 24 h降雨量與最大1h降雨量的風(fēng)險(xiǎn)聯(lián)合頻率分布圖

以宿遷市50年一遇的設(shè)計(jì)雨型為例，計(jì)算4種方法雨型的24h降雨量，最大1、3、6 h降雨量，見圖7～8，各自的風(fēng)險(xiǎn)聯(lián)合分布概率，結(jié)果見表6。在24 h設(shè)計(jì)雨量一定的條件下，時(shí)段最大降雨量越大，風(fēng)險(xiǎn)聯(lián)合分布概率越小。

由表6可知：(1)同頻率雨型和衰減指數(shù)法雨型都以P-III型曲線計(jì)算的設(shè)計(jì)降雨量為基礎(chǔ)，風(fēng)險(xiǎn)聯(lián)合分布概率相同，但同頻率雨型雨峰偏后。(2)SCS雨型雨峰較大，24 h與最大1 h對(duì)應(yīng)風(fēng)險(xiǎn)聯(lián)合分布概率為0.0001(較小)，較安全；24 h與最大6 h的風(fēng)險(xiǎn)聯(lián)合分布概率為0.010306，表明在最大6 h時(shí)段的降雨不太集中。(3)PDF法雨型雨峰與同頻率雨型一樣，但24 h和最大3、6 h的風(fēng)險(xiǎn)聯(lián)合分布概率都為0.0000001(最小)，表明最大6 h和最大3 h的降雨較集中，設(shè)計(jì)雨型偏安全。

表6 50年一遇下4種雨型(風(fēng)險(xiǎn))聯(lián)合分布概率

4 結(jié) 論

(1)暴雨衰減指數(shù)法和同頻率法原理相近，是以時(shí)段設(shè)計(jì)雨量為基礎(chǔ)來設(shè)計(jì)雨型；SCS-Ⅱ法是將美國的設(shè)計(jì)結(jié)果運(yùn)用于宿遷市來設(shè)計(jì)雨型；PDF法是利用概率密度函數(shù)來設(shè)計(jì)雨型。本文嘗試運(yùn)用該3種新方法推求宿遷市設(shè)計(jì)雨型，其雨型結(jié)果相對(duì)于同頻率雨型來說整體雨量過程相似，雨量偏差不大，結(jié)果合理，因此該方法可為宿遷市雨型的確定提供參考，也能為國內(nèi)設(shè)計(jì)雨型的研究發(fā)展提供新思路。

(2)從Copula函數(shù)中的聯(lián)合分布概率的思路出發(fā)，重新定義風(fēng)險(xiǎn)率，選用GH Copula函數(shù)計(jì)算風(fēng)險(xiǎn)聯(lián)合分布概率來進(jìn)行4種雨型的對(duì)比，來分析設(shè)計(jì)雨型的優(yōu)劣，Copula函數(shù)可以作為雨型對(duì)比分析的一種方法。