基于Copula函數的設計暴雨雨型研究

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

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

1 研究背景

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

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

2 設計雨型的研究方法

宿遷地處亞熱帶和暖溫帶的過渡地區，主要受北方西風槽、冷渦和熱帶臺風，及當地生成的江淮切變線和氣旋等天氣系統影響，暴雨主要發生在7-9月份，具有時長、量大、面廣的特點，其中特大暴雨會對其城市防洪排澇系統帶來巨大挑戰。宿遷閘站位于宿遷中心城區邊緣，離中心城區最近，是水文部門的國家基礎站，其可靠性、一致性和代表性良好。

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

2.1 同頻率法雨型

同頻率法即選取當地實測的典型暴雨，對其降雨過程進行同頻率分時段縮放[5]。放大倍比系數公式：

(1)

式中：X設計為不同重現期的設計降雨量，mm；X典型為典型暴雨降雨量，mm。

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

采用皮爾遜III型目估適線法[6]計算宿遷市各重現期各時段設計降雨量如表1。

表1 不同重現期下各歷時設計降雨量 mm

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

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

2.2 SCS雨型

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

表2 SCS-Ⅱ雨型時程分布數據

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

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

2.3 暴雨衰減指數法雨型

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

(2)

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

當Bk=1時，反推算匯流時間區間為(t1，t2)的暴雨衰減系數為：

(3)

式中：t1，t2為匯流時間，h；ht1，ht2為對應匯流時間的設計降雨量，mm。

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

表3 不同時間區間的np值

已知np的前提下，各時段設計雨量可根據以下公式計算[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)

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

圖3 宿遷市24 h設計雨型(暴雨衰減指數法)

2.4 PDF法雨型

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

(1)Gumbel分布

(7)

式中：α與β為參數。

(2)Lognormal分布

(8)

式中：η與σ為參數。

(3)Weibull分布

(9)

式中：γ與α為參數。

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

Tp=α

(10)

(11)

(12)

參數α，β要同時滿足上述的3個目標函數，顯然這是一個多元非線性方程組的求解問題，不能直接解出數值，需利用最優化算法進行求解。因此本文采用絕對值求和法，將多目標函數問題轉化成單目標函數最優化問題，只要求解出絕對值之和最小時的參數即可，參數結果應滿足下式：

(13)

2.4.2 粒子群優化算法 PDF參數求解問題轉化成了單目標函數minZ的最優化求解問題。粒子群算法(PSO)算法簡單，求解速度快，適合求解實數問題，因此本文采用PSO對函數minZ進行最優化求解，解出Z最小時對應的參數結果。

粒子群優化算法采用速度-位置模型，每個粒子代表著參數解空間的一個，解的好壞程度由適應函數的結果來決定[10]。本文的粒子群算法的學習因子取c1=c2=2。為了使粒子不遠離解空間，粒子移動速度應控制在[-vmax，vmax]，一般取vmax=k·vmax，本文取k=0.5。粒子在解空間不斷更新迭代直到達到規定的迭代次數或達到規定的誤差標準為止[11]。

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

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

(1)擬優平方和準則法

(14)

(2)擬優絕對值準則法

(15)

(3)概率點據相關系數檢驗法

(16)

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

由上述公式計算出1981-2016年優度檢驗均值，結果見表4。擬合優度檢驗中RMSE和MAE越小越優，PPCC越大越優。對于宿遷市實測數據來說，Gumbel分布是最優分布，可作為其設計雨型。

表4 擬合優度檢驗結果

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

3 基于Copula函數設計雨型的比較分析

不同設計雨型的優劣沒有固定的方法來判別，本文嘗試利用Copula函數，從設計降雨總量分別與最大1h雨量、最大3h雨量、最大6h雨量間的風險聯合分布概率入手，對上述4種雨型進行對比分析。

Copula函數是定義在[0，1]區間均勻分布的聯合概率分布函數[13]。由宿遷市36 a實測數據得出變量間為正相關，故本文選取了Gumbel-Hougaard Copula函數來分析設計雨型。GH Copula函數的表達式為：

(17)

θ∈[1,∞]

式中：θ為參數。

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

(18)

式中：(xi，yi)、(xj，yj)分別為兩組獨立同分布的觀測數據；sign為符號函數。

表5 P-III邊緣分布函數的參數值

計算得24 h降雨量與最大1、3、6 h降雨量的Kendall秩相關系數τ分別為0.5016、0.5952和0.6302，則對應GH Copula的參數θ分別為2.01、2.47和2.7。

以24 h降雨量與最大1 h降雨量的Copula概率分布關系為例來詳細介紹。采用Genest-Rivest法驗證理論聯合分布函數與經驗聯合分布函數的擬合程度，結果如圖6所示。各點均勻的分布在45°線左右，則GH Copula具有良好的擬合度[15-16]。

本文定義風險率為滿足24 h設計降雨量條件下，最大1h(或3、 6 h)降雨量超出對應設計時段降雨量發生的條件概率[5]，計算公式如下式：

(1)

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

圖4 2000年3個分布函數的擬合曲線圖

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

圖5 宿遷市24 h設計雨型(PDF法雨型)

圖6 經驗聯合分布函數與理論聯合分布函數

圖7 24 h降雨量與最大1h降雨量的聯合頻率分布圖

圖8 24 h降雨量與最大1h降雨量的風險聯合頻率分布圖

以宿遷市50年一遇的設計雨型為例，計算4種方法雨型的24h降雨量，最大1、3、6 h降雨量，見圖7～8，各自的風險聯合分布概率，結果見表6。在24 h設計雨量一定的條件下，時段最大降雨量越大，風險聯合分布概率越小。

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

表6 50年一遇下4種雨型(風險)聯合分布概率

4 結 論

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

(2)從Copula函數中的聯合分布概率的思路出發，重新定義風險率，選用GH Copula函數計算風險聯合分布概率來進行4種雨型的對比，來分析設計雨型的優劣，Copula函數可以作為雨型對比分析的一種方法。