基于Kendall重現期的華南中小流域洪水峰量聯合分布研究

趙玲玲, 楊 興, 劉麗紅, 劉昌明

(1.廣東省遙感與地理信息系統應用重點實驗室, 廣東 廣州 510070;2.廣州地理研究所 廣東省地理空間信息技術與應用公共實驗室， 廣東 廣州 510070；3.廣州地理研究所， 廣東 廣州 510070； 4.中國科學院 地理科學與資源研究所， 北京 100101；5.安徽理工大學 地球與環境學院， 安徽 淮南 232000； 6.北京師范大學 水科學學院， 北京 100875)

水文極端事件的風險概率及其相應的重現期標準是當前設計洪水計算中迫切需要解決的重要科學與工程應用問題。受亞熱帶季風性氣候影響，華南地區天氣系統復雜中小尺度天氣系統十分活躍，汛期降雨具有強度大，季節性強，時間短、范圍小、致災性重等特點[1]，山區地貌陡峻，河道坡降大，山洪災害頻發。在氣候變化背景下，山區中小流域增大了極端降水頻率與強度，產生的山洪災害已成為防災減災的重點[2]。研究年內洪水最大洪峰和洪量的分布類型和重現期，均可為防洪工程的建設和優化提供依據。但是單變量頻率的分析并不能滿足洪水多個特征屬性的特點，且變量之間存在相依性，因此同時考慮洪峰和洪量兩個洪水特征之間的聯合頻率，對流域洪水特征的分析是一個重要的研究方面。近年來，對洪峰、洪量和洪水歷時等多變量聯合分析受到關注[3-11]，其中Copula函數能夠靈活的聯結兩變量或多變量的水文因子聯合分布，在水文變量分析中得到廣泛的應用。Nelson[12]總結了Copula相關研究領域的主要成果，謝華等[13]對水文頻率研究中二維聯合概率采用的4種阿基米德族Copula函數及其綜合選優做了述評。侯蕓蕓等[14]、陳子燊等[9]應用Copula函數，探討洪水三變量的聯合概率分布和條件概率分布，驗證了Copula函數可以在三維洪水變量中的適用。郭生練等[15-16]對多變量水文分析計算中的應用與研究進展作了述評，指出對于給定的重現期如何合理地選擇聯合設計值是關鍵問題。

目前聯合頻率分析多集中在對兩要素的“或”和“且”重現期上，Salvadori等[17]和有關研究[9,15]表明使用Copula函數計算多變量“或”聯合重現期和“且”聯合重現期十分簡便，可為風險分析提供一種非常簡單而又有效的方法。但“或”聯合重現期和“且”聯合重現期在危險域或安全域劃分上存在局限性，Salvadori等[18]等引入了一個新的可與特定事件聯合重現期相關聯的分布函數——Kendall測度，并定義了一個新重現期Kendall重現期，其涵義為超過閾值事件的平均到達時間(臨界事件)。相比于傳統重現期，Kendall聯合重現期的提出改進了多變量聯合設計的可靠性[19]，為處理潛在危險(破壞性)的隨機事件的頻率分析領域提供了新的研究途徑，Corbella S等[20]等在海岸侵蝕中進行了應用；陳子燊等[8]在城市洪澇中基于Kendall重現期推算的不同歷時暴雨組合的設計暴雨分位值，驗證了Kendall重現期優于傳統的重現期；劉章君等[21]考慮了洪峰、洪量與水庫調洪規則的交互作用中發現Kendall和生存Kendall重現期也存在不同程度的偏低或偏高；范嘉煒等[22]、史黎翔等[23]在城市洪澇、無定河流域中應用研究，也驗證了Kendall重現期的合理性。但是對于山區中小流域山洪特點下Kendall重現期的相關研究存在不足，能否適用于山區中小流域洪水設計工程中還有待研究。因此，本文擬通過對華南三個山區中小流域的洪峰和洪量聯合分布的實例，分析“或”重現期、“且”重現期和Kendall重現期的設計水平之間的差異，以期為多變量洪水頻率分析在山區洪水中的應用提供理論參考。

1 理論與方法

1.1 Copula函數

Copula理論指出多個變量的聯合分布可分解為多個不同的邊緣分布，通過一個Copula函數構建聯合分布。據此理論，設(x,y)為二維隨機變量，u=FX(x)，v=FY(y)為連續的邊緣分布函數，則有唯一的Copula函數C使得：

F(x,y)=P(X≤x，Y≤y)

=C〔FX(x)，FY(y)〕=(u,v)

(1)

式中：F(x,y)為聯合分布函數，關于Copula函數性質的更詳細描述可參閱Salvadori等(2007)[24]的論述。

1.2 重現期定義

(2)

(3)

式中：μ——兩個連續事件的平均到達時間。

1.2.2 Kendall分布函數與Kendall重現期 從“或”和“且”重現期可知，相同的累積頻率均可產生相同的重現期，不會因u，v組合事件的不同而改變。為解決由“或”和“且”重現期存在的對安全事件與危險事件錯誤的識別，又將其分為亞臨界(安全域)、臨界(警戒事件)和超臨界(危險域)3種情景[24-25]。通過判定累積概率不超過某臨界概率，將多維的極值事件投射為一維分布，完全區分了安全事件與危險事件在空間域的分布。黃強等[11]對此作了詳細的圖解說明。與Copula函數累積概率為t的(u,v)組合值相關聯的Kendall測度KC為：

(4)

式中：φ(t)為KC生成元；φ′(t)為φ(t)的右導數。由Kendall測度確定的重現期稱為Kendall重現期(TKEN)：

(5)

Graler等[26]對二維Copulas的“或”、“且”和Kendall重現期三種聯合重現期(JRP)的不等關系作了圖示解釋(圖1)：對于一固定的設計事件(u,v)，其累積分布單位平方圖內不同的重現期TOR，TKEN和TAND可以用1/(1-面積(安全事件))表示。如圖1所示，“或”重現期定義僅將左下角矩形中的所有事件視為安全的。Kendall重現期將左上角和右下角的曲線區域(KEN)與左下角的矩形區劃為安全域，從而使同一設計事件(u,v)產生的重現期比“或”重現期更大。“且”重現期則進一步添加了左上角和右下角的矩形，從而得到最大的重現期。由式(2)—(5)和圖1可知“或”、“且”和Kendall重現期之間的不等關系為：TOR≤TKEN≤TAND。

注：(u,v)為某一設計事件；C(u,v)為Copula水平曲線；KEN為Kendall重現期定義的安全域；AND為And重現期定義的安全域；OR為or重現期定義的安全域。

圖1 3種聯合分布重現期定義的圖示說明[26]

1.3 水文特征設計值

由概率論，Q≥q條件下W≥w的累積概率分布為:

P(W≥w|Q≥q)=

(6)

此條件概率屬于超值概率，可定義為遭遇概率。分析二者的遭遇概率有助于進一步認識流域洪水事件的統計規律，對于防洪管理決策具有重要的參考意義。

由于對某預定的重現期存在無數個滿足防洪標準的多變量分位值組合，如何合理地推算聯合設計值的問題成為一個關鍵問題[15]。有關研究指出，在重現期分位值相同情況下，必然存在一個組合值使得聯合概率密度達到最大值，可利用以下公式推算聯合分布設計值：

(7)

f(u,v)=c(u,v)f(u)f(v)

(8)

式中：(um,vm)為兩變量聯合概率密度 達到最大值時的組合設計值；f(u)，f(v)為邊緣分布的概率密度函數；c(u,v)為二維Archimedean Copula的概率密度函數。

2 實例分析

2.1 流域概況

本文選取廣東省中小流域曹江、田頭水和羅壩水流域作為研究區，各站點位置為大拜水文站E111°09′N22°03′，結龍灣水文站E114°11′N24°54′，赤溪水文站E113°08′N25°03′，其中田頭水流域只有廣東省內水系。3個流域均為山區中小流域，地形崎嶇、河流坡度大并且人為影響較小，其中曹江流域是廣東省三大暴雨中心，出口斷面大拜水文站集水面積394 km2，流域多年平均年雨量可達2 160 mm，最大平均年雨量達3 150 mm。赤溪水文站位于廣東樂昌市慶云鎮赤溪，是珠江流域北江水系武江支流田頭水流域的控制站，集水面積為396 km2。結龍灣水文站是珠江流域北江水系一級支流墨江下游羅壩水流域的控制站，集水面積81 km2，流域多年平均年徑流量2.75×108m3。

本研究分別收集了大拜水文站1967—2013年、赤溪水文站1967—2016年、結龍灣水文站1967—2016年場次洪峰流量數據。依據超定量法[27]，首先根據3個水文站枯水年最大洪峰流量為最小閾值，然后按各場洪水流量過程線提取洪峰流量并計算出該場洪水的洪量與歷時，并剔除較差線型的洪水樣本。為保證不同洪水之間的獨立性，各場次洪峰發生的時間間隔要求大于流域匯流時間。3個流域分別抽取了231，169，264場洪峰流量和相應洪量作為峰量聯合分布的分析樣本。

根據洪水流量過程線(圖2)提取洪峰流量Q(m3/s)、洪水歷時D(天)和洪水總量(洪量)W：

D=te-ts

(9)

(10)

式中：ts為洪水開始時間；te為洪水結束時間；q為日流量序列。

注：S，e為洪水基流切割點；qte為洪水開始點流量；qts為洪水結束點流量；Q為流量；D為洪水歷時。

圖2 洪水過程及相應的洪量與歷時

2.2 Copula函數的選取

根據3個水文站多年洪水流量數據，用水文頻率分析中常用的4種三參數概率分布：對數正態分布(GNO)、廣義極值分布(GEV)、廣義Logist(GLO)分布和皮爾遜三型分布(P-Ⅲ)，使用較穩健的線性矩方法[28]，分別對洪峰流量和洪水總量樣本加以擬合。并對4種參數擬合的結果采用均方根誤差(RMSE)和概率點據相關系數(PPCC)檢驗其擬合優度選擇最優水文頻率分布函數，擬合結果見表1。大拜站洪峰洪量序列宜選用GNO分布，赤溪站分別選用GLO分布和GEV分布，結龍灣站分別選用GLO分布和P-Ⅲ分布。

表1 3個流域洪峰和洪量的概率分布參數與擬合優度檢驗值

注：GNO為對數正態分布；GEV為廣義極值分布；GLO為廣義Logist分布；P-Ⅲ為皮爾遜三型分布；RMSE為均方根誤差；PPCC為概率點據相關系數。

計算洪峰流量和洪水總量之間的Kendall相關系數，得到大拜水文站、赤溪水文站和結龍灣水文站的洪峰Q和洪量W的Kendall相關系數分別為0.79，0.80，0.76，表明洪峰流量和洪水總量之間具有較強的相關性。

采用基于秩相關的Kendall相關系數的計算洪峰Q和洪量W聯合分布的4種阿基米德Copula參數θ，并采用Akaike信息準則法(AIC)和最小二乘法準則(OLS)確定結果(見表2)，可見4種Copula函數擬合結果較好，依據AIC和OLS結果選擇最優Gumbel Copula函數來構建3個流域洪峰Q和洪量W的聯合分布函數，各站構建的Copula分布模式如下：

大拜站：C〔FQ(q),FW(w))〕=exp{-〔(1-lnFQ(q))4.703+(-lnFW(w))4.703〕1/4.703

(11)

赤溪站：C〔FQ(q),FW(w))〕=exp{-〔(1-lnFQ(q))4.967+(-lnFW(w))4.967〕1/4.967

(12)

結龍灣站：C〔FQ(q),FW(w))〕=exp{-〔(1-lnFQ(q))4.142+(-lnFW(w))4.142〕1/4.142

(13)

表2 4個Copula函數的參數及其擬合優度指標

注：OLS為最小二乘法準則，AIC為信息準則法，θ為參數。

2.3 條件概率分布

洪峰流量和洪水總量都屬于洪水過程的隨機變量，在最大洪峰流量Q和最大洪水總量W都屬于洪水過程的隨機變量，估計出洪峰Q和洪量W之間的聯合分布函數之后，就可以推求出給定某種事件發生概率的條件下的另一種事件發生的概率。分析特定設計洪峰流量條件下出現洪水總量的概率分布，由公式(6)可以求出洪峰洪量條件概率值。

根據洪峰洪量條件概率結果(見表3)，3個流域出現概率大于等于表中概率的洪峰流量時，洪水峰量同頻率遭遇概率很大，且由于洪峰洪量具有高相關性，還有很大可能出現洪峰遭遇更小頻率的洪量。其中，曹江流域大拜站峰量同頻率遭遇大于84.1%，主對角線以上二者遭遇的概率則大于98.4%；田頭水流域赤溪站峰量同頻率遭遇大于85%，主對角線以上二者遭遇的概率則大于98.7%；羅壩水流域結龍灣站峰量同頻率遭遇大于81.8%，主對角線以上二者遭遇的概率則大于97.3%。

由以上分析可見，此3個流域洪水峰量的遭遇風險概率基本接近，當洪峰大于或等于某一設定頻率時，洪量出現大于該頻率的條件概率隨之增大。上述結果表明洪水峰量聯合分布可能存在著多種結果，分析多種組合出現的不同遭遇概率有利于現實中防汛減災的風險管理。

表3 華南中小流域洪峰洪量條件概率〔P(W≥w|Q≥q)〕

2.4 聯合分布重現期

采用Gumbel Copula聯結函數根據式(7)求洪峰Q和洪量W組合的“或”、“且”和二次重現期結果(見表4)。從3個流域洪水峰量聯合分布計算結果可看出，對于特定的設計重現期，Kendall重現期都介于“或”和“且”重現期之間，并且大于設定重現期。“或”聯合重現期最小且小于設定的重現期，“且”重現期最大且大于設定的重現期。“或”、“且”和二次重現期聯合分布重現期與特定的設計重現期之間兩者的相對為誤差為，TOR：-0.12%～-0.15%，TAND：0.16%～0.22%，TKEN：0.08%～0.11%。Kendall重現期聯合分布設計值最為接近設定重現期，對洪水峰量之間任一要素超標致災的重現期標準宜采用Kendall重現期。3種聯合重現期較為接近，或與此3個中小流域同處華南氣候區，洪水事件具有相近的背景有關。

表4 華南中小流域洪水峰量聯合分布的重現期

2.5 設計洪水分位值

根據重現期計算方法的不同，分別求設定重現期下各水文站的洪峰流量和洪水總量單變量設計值列于表5。分析表5中3個流域不同重現期下洪峰流量與洪水用量的設計值特征，可以得出以下結論。

(1) 對于設定的重現期，3個流域的單變量設計洪峰流量和設計洪水總量及其聯合設計值，以赤溪站最大，大拜站次之，結龍灣站最小。雖然赤溪站所在的田頭水流域和大拜站所在的曹江流域集水面積大致相同，但主要與田頭水流域2006年7月14—15日2 d內遭遇652.5 mm特大降雨形成的超百年一遇大洪水有關。

表5 3個流域不同重現期下洪峰流量與洪水總量的設計值

注：Q為洪峰流量(m3/s)；W為洪水總量(106m3)。

(2) Kendall重現期推算的值小于“或”重現期設計值和邊緣分布推算的設計洪峰流量和設計洪水總量設計值，兩者的相對為誤差為，大拜站：洪峰流量:-0.7%～1.1%，洪量:-1.3%～0.7%；赤溪站：洪峰流量:-4%～0.4%，洪量:0%～3%；結龍灣站：洪峰流量:-5.9%～-3.4%，洪量:-2.2%～-0.1%。

上述結果表明，Kendall重現期推求的設計洪水值大于洪水、峰量邊緣分布設計值，并介于“或”和“且”重現期設計值之間。與表4聯合分布結果一致，因此Kendall重現期可以作為洪水峰量聯合分布的待選方法，為防洪工程安全與風險管理提供更好的選擇。為進一步比較分析，推算洪峰洪量同頻率分布設計值[29]：

u1=u2=〔1-(1-Tu1,u2)〕α；

Q=F-1(u1)；W=F-1(u2)

(14)

式中：α=2-1/θ；Tu1,u2為“或”重現期；F-1(ui)為邊緣分布函數的反函數。

3個流域不同重現期下洪峰流量與洪水總量的設計值(見表5)，進一步比較分析洪水峰量聯合分布的同頻率設計值和聯合概率密度最大值推算的重現期設計值，可以看出設定重現期下兩者設計值十分接近，并且均低于“或”重現期設計值。表明了對于工程設計來說“或”重現期或同頻率分布設計值存在“高估”的問題，會提高水工建筑建設成本。

一場洪水的造成的損失并不是單個特征變量決定的，而是多個洪水特征變量共同作用的結果，例如同一量級的洪水并不一定造成同樣的影響，不同大小的洪量也可能導致同樣程度的破壞，因此多變量聯合分布的研究十分重要。其中一個重要方面是變量之間的建模，Copula聯結函數在多元洪水特征分析方面的潛力；另一個重要方面是重現期方法的選取，不同方法的選擇直接影響到建造結構的安全和成本。相對于傳統重現期，Kendall重現期設計值介于“或”與“且”重現期之間，且略高于設計重現期。Kendall重現期在工程安全設計理念高于“或”重現期設計值，對于工程經濟理念來說高于“且”重現期設計值。上述洪峰和洪量的Kendall相關性較高，且主對角線上條件概率超過81%表明洪水峰量聯合分布存在多種結果。由表4和表5分析可知，Kendall重現期在3個流域有較好的應用結果。且黃強等[11]研究也表明傳統的“或”和“且”多變量重現期對安全與危險域的識別具有局限性。因此，可以說明，Kendall重現期在洪水峰量聯合分布計算中優于傳統重現期，可以為山區中小流域防洪工程安全提供更多的選擇。

3 結 論

本文對比分析了廣東省中小流域曹江、田頭水和羅壩水流域洪峰洪量之間的聯合分布及其重現水平，有以下結論：

(1) 由4種累積分布函數擇優構成了3個流域不同的洪水峰量聯合，AIC和OLS擬合結果表明使用Gumbel Copula構建了最佳的洪峰洪量聯合分布;洪峰Q和洪量W的Kendall相關系數分別為0.79，0.80，0.76,表明洪峰流量和洪水總量之間具有較強的相關性。

(2) 作為山區暴雨洪水成因的洪水過程,3個流域洪水峰量相關性高,主對角線以上的條件概率超過81%,洪水峰量遭遇風險概率大且基本接近,分析多種洪水峰量組合出現的不同遭遇概率有利于防汛減災的風險管理。

(3) 相對于“或”重現期,采用Kendall測度計算的Kendall重現期可更好地區分超臨界事件的風險率。Kendall重現期推算的洪峰洪量設計洪水值介于“或”重現期與“且”重現期設計值之間,接近于邊緣分布設計值。在經濟安全方面,Kendall重現期設計洪水更優于傳統重現期設計方法,可為防洪工程風險管理與設計提供新的選擇與參考依據。