蔡奇霖，徐丹麗，王兆禮，賴成光*

(1.華南理工大學土木與交通學院，廣東 廣州 510640；2.廣東省水利電力勘測設計研究院有限公司，廣東 廣州 510635)

在全球氣候變暖的背景下，不少地區(qū)的臺風與暴雨頻度、規(guī)模、時空特征等均發(fā)生了較大變化[1-3]，這些變化給地區(qū)水庫設計、建設、管理等帶來諸多挑戰(zhàn)。在水庫大壩設計過程中，確定合理的設計風速對雍水及風浪爬高計算至關重要，決定了大壩的綜合超高以及安全裕度。然而，當前中國大部分技術規(guī)范僅考慮單一的風速變量確定設計風速，鮮有定量考慮強風和暴雨遭遇的情況，如SL 274—2020《碾壓式土石壩設計規(guī)范》和SL 25—2006《砌石壩設計規(guī)范》是按大壩的重要性給出相應的風速標準，沒有考慮風速與暴雨的遭遇。盡管SL189—2013《小型水利水電工程碾壓式土石壩設計規(guī)范》要求在沿海地區(qū)應該考慮最大風力與暴雨同時出現(xiàn)的條件，但如何執(zhí)行也沒有明確規(guī)定。廣東省是中國受臺風和暴雨影響最嚴重的沿海地區(qū)之一，雙臺風甚至是三臺風同時出現(xiàn)的現(xiàn)象也并非歷史罕見。已有研究表明該地區(qū)登陸臺風呈現(xiàn)出頻數(shù)減少但強度明顯增強的特點[4-6]，而降雨則表現(xiàn)出強降水增加，弱降水減少的特點[7-9]。在降雨和風速均發(fā)生改變的背景下，當二者發(fā)生遭遇時，按照目前有關規(guī)范計算的壩頂超高(特別是沿海地區(qū)的水庫)是否仍能滿足安全要求，這是一個值得深入研究的課題。

在眾多的多變量分析方法中，Copula函數(shù)因其不受邊緣分布限制的優(yōu)點，已在水文頻率分析中得到了廣泛應用[10-11]。徐明等[12]利用Copula函數(shù)針對上海地區(qū)建立了日降雨量極值與極大風速極值的風雨聯(lián)合分布模型，結(jié)果表明聯(lián)合分布模型得到的風雨同現(xiàn)超越概率能較好地反映風險程度。侯靜惟等[13]以海南省計值風速和總降水量為研究樣本，建立風雨聯(lián)合分布模型，并分析了“或”和“且”2種重現(xiàn)期的風雨重現(xiàn)水平。Dong等[14]針對崇明島地區(qū)基于Copula函數(shù)建立了風速和雨強的聯(lián)合分布模型，并通過了實測數(shù)據(jù)檢驗。徐宗學等[15]基于Copula函數(shù)對深圳市洪潮組合進行風險分析，結(jié)果表明降水與潮位的雙閾值組合風險率小于單閾值組合風險率，設計時應更加關注降水和潮位的雙閾值組合風險率。馬川惠等[16]基于Copula函數(shù)針對涇河流域建立水沙聯(lián)合分布模型得到了兩變量聯(lián)合設計值的最可能模式。李天元等[17]以隔河巖水庫為例，基于Copula函數(shù)構建洪峰、洪量間的二位聯(lián)合分布，推導了兩變量同頻率組合和條件期望組合兩種聯(lián)合設計值組合。林嫻等[18]以武江流域?qū)崪y資料為基礎，基于Copula函數(shù)構建了洪峰、洪量與洪水歷時間的聯(lián)合概率分布，結(jié)果表明同頻率下聯(lián)合分布的洪水設計值相對于單變量設計值偏安全。曾明等[19]基于上海徐家匯站的逐日降雨資料，采用Copula函數(shù)構建了年最大1 d與3 d降水量的聯(lián)合分布模型，結(jié)果表明考慮最大1 d和3 d降水量的遭遇組合，有利于提升防洪排澇能力。綜上可知，Copula函數(shù)已被廣泛用來構建聯(lián)合分布連接函數(shù)，能夠有效地描述水文事件的內(nèi)在規(guī)律和特征屬性之間的相互關系。

盡管目前關于風速和降雨的聯(lián)合分布研究已有較多成果，但關于廣東沿海地區(qū)大型水庫風速與降雨的聯(lián)合分布研究尚不多見?；诖耍疚囊詮V東省臺山大隆洞水庫為研究案例，基于Archimedean Copula函數(shù)和Kendall分布函數(shù)分別構建了年最大1 d降雨量、年最大3 d降雨量、年最大7 d降雨量及其對應時段最大風速的二維聯(lián)合分布函數(shù)，推求不同時段最大降雨量以及最大風速在“或”“且”和Kendall重現(xiàn)期的重現(xiàn)期水平，并對比分析3種遭遇情景結(jié)果間的差異，最后基于遭遇分析的風速結(jié)果與當前技術規(guī)范進行對比，并提出相關建議。研究結(jié)果以期為沿海地區(qū)水庫大壩的設計、建設以及管理提供科學的參考依據(jù)。

1 水庫簡介

大隆洞水庫位于廣東省臺山市東南部，是廣東省受臺風、暴雨影響最為頻繁的大型水庫之一。由大隆洞水庫大壩基本資料登記表(BHF24000112)和水庫管理處有關資料獲知，該水庫于1959年建成，2003年經(jīng)加固改造后水庫集雨面積達148 km2，水庫總庫容2.96億m3，為大(2)型水庫，設計重現(xiàn)期為100年一遇。水庫主壩壩型為均質(zhì)土壩，壩頂高程(不含防浪墻)39.2 m，壩頂長度為441 m，壩頂寬度6 m，壩基采用混凝土防滲墻防滲；水庫正常蓄水位30.80 m，設計洪水位35.42 m，校核洪水位38.11 m。

降雨和風速數(shù)據(jù)來源于水庫附近的臺山站(112°78′E、22°25′N)，其中降雨(逐日)數(shù)據(jù)序列長度為1960—2018年，最大風速(逐日)數(shù)據(jù)序列長度為1973—2018年。考慮到規(guī)范要求和歷史設計資料一致性問題，采用年最大值法(AM)提取每年的最大日降雨量和最大風速作為年最大日降雨量和年最大風速，并重新組成年最大日降雨量和年最大風速序列，二者其變化趨勢見圖1。此外，分別使用滑動法獲得的年最大3 d和7 d降雨量并形成新的暴雨序列，并取對應時間段內(nèi)日最大風速作為風速序列，以分析二者之間的遭遇規(guī)律。

a)年最大日降雨

b)年最大風速

2 研究方法

首先分別采用P-III、Gamma和GEV 3種邊緣分布函數(shù)對3組降雨和1組風速序列進行擬合，并采用K-S檢驗和PDE方法選出最優(yōu)邊緣分布函數(shù)。然后采用Gumbel、Frank和Clayton Copula 3種Copula函數(shù)構建強風暴雨的聯(lián)合概率模型，并采用OLS和AIC準則選出最優(yōu)擬合函數(shù)。最后分別計算風雨兩變量的“或”“且”和Kendall重現(xiàn)期設計值，并基于風速分析結(jié)果與現(xiàn)行技術規(guī)范進行對比分析。

2.1 邊緣分布函數(shù)

邊緣分布函數(shù)均已在風速、暴雨、風暴潮、洪水等序列的極值分布規(guī)律分析方面取得了理想的效果[20]。通過綜合對比和考慮，最終選擇皮爾遜III型(P-III)、伽馬分布(Gamma)、廣義極值分布(GEV)3種常用函數(shù)，各函數(shù)見表1。

表1 獨立概率分布函數(shù)與概率密度函數(shù)

邊緣分布采用概率分布誤差率(Probability Distribution Error，PDE)[21]和Kolmogorov-Smirnov(K-S)檢驗[22]進行擬合優(yōu)度評價，以確定最優(yōu)邊緣分布函數(shù)。其中概率分布誤差率計算公式如下：

(1)

式中Fn(x)——經(jīng)驗分布函數(shù)；F(x)——理論分布函數(shù)；n——樣本的個數(shù)。

2.2 Copula函數(shù)

水文分析中常用的Copula函數(shù)包括經(jīng)驗Copula函數(shù)、橢圓Copula函數(shù)及Archimedean Copula函數(shù)等，其中Archimedean Copula函數(shù)因其計算簡便、形式多樣，已在實際應用中得到了廣泛應用，效果良好[23-24]。本研究將選用Gumbel、Frank和Clayton Copula 3種函數(shù)(表2)，聯(lián)合分布采用AIC(Akaike information criterion)信息準則法[25]及OLS(Ordinary least squares)離差平方和最小準則[26]進行擬合優(yōu)度評價，AIC值與OLS值越小則表明擬合效果越好。

表2 常見Archimedean Copula函數(shù)

OLS準則是通過計算理論值和實測值的均方根誤差RMSE來定量地評估擬合誤差，其計算公式為：

(2)

式中Femp(xi1,xi2…,xin)——經(jīng)驗頻率值；C(ui1,ui2…,uin)——理論頻率值。

AIC指標用于對Copula函數(shù)擬合情況進行評價，其計算公式為：

C(ui1,ui2…,uin)]2}+2k

(3)

式中n——函數(shù)維數(shù)；k——模型參數(shù)個數(shù)。

2.3 重現(xiàn)期計算

水文分析中常用的重現(xiàn)期包括“或”重現(xiàn)期和“且”重現(xiàn)期，“或”重現(xiàn)期表示2個變量中至少有1個變量超過特定頻率時的重現(xiàn)期，“且”重現(xiàn)期表示2個變量同時超過特定頻率時的重現(xiàn)期[27]?！盎颉敝噩F(xiàn)期和“且”重現(xiàn)期分別用式(4)、(5)表示：

(4)

(5)

式中P——同現(xiàn)、聯(lián)合概率；X1、Y1——變量序列；Fx1、Fx2——邊緣分布函數(shù)；F——聯(lián)合分布函數(shù)；Tor、Tand——“或”“且”重現(xiàn)期。

Kendall重現(xiàn)期考慮了曲線概率等值情況，相比于“或”“且”重現(xiàn)期更為嚴謹[28]。該重現(xiàn)期計算中引入了Kendall測度KC來計算，其中KC與Copula函數(shù)有關，計算公式如下：

(6)

(7)

3 結(jié)果分析

3.1 相關性與遭遇統(tǒng)計

分別計算1973—2018年最大風速與最大日降雨量、最大3 d降雨量、最大7 d降雨量序列的相關性，可得其相關系數(shù)分別為0.153、0.012和-0.036，表明3組序列相關性極弱。然而，相關性分析難以獲悉二者的實際遭遇情況。分別選取年最大風速以及最大日降雨、最大3 d和最大7 d降雨統(tǒng)計大隆洞水庫的強風暴雨遭遇情況(表3)。在46年的樣本中，當大隆洞水庫發(fā)生年最大風速時剛好遭遇最大日降雨的次數(shù)為6次，占比為13.04%；當發(fā)生年最大風速剛好遭遇最大3天降雨的次數(shù)為8次，占比為17.39%；當發(fā)生年最大風速剛好遭遇最大7 d降雨的次數(shù)為7次，占比為15.22%。以上表明，盡管大隆洞水庫的風速和暴雨序列相關性極弱，但出現(xiàn)年最大風速和年最大暴雨遭遇的情況并非罕見。

表3 大隆洞水庫年最大風速與不同時段最大降雨量遭遇情況

3.2 邊緣分布及Copula函數(shù)擬合優(yōu)度評價

以1973—2018年間年最大日降雨、最大3 d降雨、最大7 d降雨以及對應時段內(nèi)的日最大風速為樣本，采用P-III、GEV和Gamma 3種函數(shù)進行邊緣分布擬合，邊緣分布采用K-S、PDE檢驗進行擬合優(yōu)度評價(表4)。

表4 各序列邊緣分布的參數(shù)估計及擬合優(yōu)度

結(jié)果表明，最大日降雨和對應日最大風速的最優(yōu)函數(shù)均為GEV分布，最大3 d降雨和對應日最大風速的最優(yōu)函數(shù)分別為P-III分布和GEV分布，最大7 d降雨和對應日最大風速的最優(yōu)函數(shù)均為GEV分布。選擇Gumbel、Frank及Clayton Copula 3種函數(shù)分別構建不同序列間的聯(lián)合分布模型，采用AIC信息準則法和OLS準則進行擬合優(yōu)度評價，Copula函數(shù)擬合優(yōu)度與OLS、AIC值呈負相關關系，結(jié)果見表5。

表5 不同序列Copula聯(lián)合分布函數(shù)擬合優(yōu)度

結(jié)果表明，日最大風速與最大日降雨、最大3 d降雨、最大7 d降雨聯(lián)合分布的最優(yōu)函數(shù)均為Clayton Copula函數(shù)，圖2為不同序列間聯(lián)合概率分布圖。隨著降雨量上限從最大日降雨的300 mm增加至最大3 d降雨的500 mm時，風速上限有一定增大，但增幅較??；隨著降雨量上限繼續(xù)增加至最大7 d降雨的600 mm時，風速上限反而呈現(xiàn)一定減小，表明風速最大值與降雨量最大值之間并沒有顯著正相關關系。

a)風速與最大日降雨

3.3 風速和降雨量設計值

3種遭遇情景下，采用等概率法分別計算200、100、50、20、10、5年一遇條件下的最大風速和最大降雨量單變量設計值，“或”“且”與Kendall重現(xiàn)期的結(jié)果見表6。結(jié)果表明，最大風速與最大日降雨量遭遇情景下，與單變量的邊緣分布結(jié)果相比，按兩變量“或”重現(xiàn)期推算的最大風速與降雨量設計值的相對差異范圍分別為25.1%～25.7%、5.8%～12.7%，高于邊緣分布設計值。按兩變量“且”重現(xiàn)期推算的最大風速與降雨量設計值的相對差異范圍分別為26.1%～56.6%、17.6%～24.4%，低于邊緣分布設計值。而按兩變量Kendall重現(xiàn)期推算的最大風速與降雨量設計值的相對差異范圍分別為18.9%～32.1%、10.7%～12.1%，低于邊緣分布設計值?？梢?，相比于單變量風速與降雨量設計值，按風速和降雨量聯(lián)合分布的“且”、Kendall重現(xiàn)期推算的風速和降雨量低于單變量的邊緣分布結(jié)果，而采用“或”重現(xiàn)期推算的風速和降雨量則高于單變量的邊緣分布結(jié)果。最大3、7 d降雨遭遇條件下具有相似的結(jié)論。

表6 不同取樣條件下不同重現(xiàn)期設計值

3.4 與當前技術規(guī)范的差異對比

前已述及，當前大多數(shù)技術規(guī)范沒有定量考慮強風和暴雨的遭遇問題，可能會對水庫的安全性產(chǎn)生一定的威脅。例如，SL 274—2020《碾壓式土石壩設計規(guī)范》中壩頂超高的計算所采用的設計風速是多年平均最大風速的1.5～2.0倍。根據(jù)風速資料計算可知，大隆洞水庫在1973—2018年的多年平均最大風速為13.91 m/s，其2.0倍則為27.83 m/s。然而，根據(jù)表6可知，最大風速與最大日降雨量遭遇情景下，等概率發(fā)生條件下100年一遇的“或”聯(lián)合分布最大風速值為36.1 m/s，顯然比2.0倍多年平均最大風速還要大29.7%。根據(jù)SL 274—2020《碾壓式土石壩設計規(guī)范》有關壩頂高程的計算，當取2.0倍多年平均最大風速時，大隆洞水庫主壩高程(壩頂+防浪墻)為39.85 m，而該水庫主壩實際高程則為40.20 m，表明按照當前的規(guī)范計算能滿足安全要求；然而，當采用100年一遇的“或”聯(lián)合分布最大風速36.1 m/s時，其主壩高程則為40.71 m，顯然超出了實際高程(表7)。

表7 計算壩頂高程與實際高程對比

4 討論

由上述分析可知，當前大隆洞水庫大壩高程是滿足現(xiàn)有規(guī)范要求的。根據(jù)調(diào)查得知，大隆洞水庫在最近幾十年汛期和臺風期間未曾出現(xiàn)水位過高導致的險情，這是由于除了高程滿足要求外，還得益于近年來天氣預報準確性不斷提高以及水庫的科學預警和調(diào)度管理，使水庫在臺風來臨前降至安全水位，大大降低了波浪爬高和水位雍高。然而，由于大隆洞水庫位于廣東沿海地區(qū)，由前述可知該地區(qū)的登陸臺風呈現(xiàn)頻數(shù)減少但強度明顯增強的特點，而降雨則表現(xiàn)出強降水增加，弱降水減少的特點。圖1也表明該水庫年最大降雨極值呈顯著上升趨勢(p<0.05)，而盡管最大風速整體呈下降趨勢，但在2000年以后也呈上升趨勢。在未來氣候變化背景下，強風和暴雨的強度可能會進一步增加，二者遭遇的概率有可能會提高。當前的技術規(guī)范僅考慮風速這一單變量，在沒有考慮二者遭遇情況下，僅采取2.0倍多年平均最大風速所計算出來的壩頂高程有可能難以滿足安全要求。因此，考慮到未來氣候變化，建議有關技術規(guī)范可考慮增加強風和暴雨遭遇的情況，適當考慮“或”聯(lián)合分布的分析結(jié)果；或者設計人員在應用當前相關規(guī)范的基礎上適當增加壩頂高程的安全余度，盡量滿足二者遭遇時的高程要求，確保水庫大壩的安全。

5 結(jié)論

以大隆洞水庫為研究實例，采用Copula函數(shù)構建了強風-暴雨聯(lián)合概率分布模型，推算了不同遭遇條件不同重現(xiàn)期的設計風速與降雨量設計值，并基于遭遇分析的風速結(jié)果與當前技術規(guī)范進行對比，得出如下主要結(jié)論。

a)大隆洞水庫年最大日降雨呈顯著上升趨勢(p<0.05)，年最大風速則呈總體下降趨勢，但2000年以后呈上升趨勢；盡管水庫的風速和暴雨序列的相關性極弱，但曾出現(xiàn)年最大風速和年最大暴雨遭遇的情況。

b)相比于單變量風速與降雨量設計值，基于Copula函數(shù)推算的結(jié)果中，聯(lián)合分布的“且”和Kendall重現(xiàn)期推算的風速和降雨量低于單變量邊緣分布的結(jié)果，而采用“或”重現(xiàn)期推算的風速和降雨量則高于單變量邊緣分布的結(jié)果。

c)按照當前規(guī)范要求，當取2.0倍多年平均最大風速時所計算出來的高程低于大隆洞水庫主壩實際高程，但當采用100年一遇的“或”聯(lián)合分布的最大風速所計算出來的主壩高程比實際高程要高；考慮到未來氣候變化，有關技術規(guī)范可考慮增加強風和暴雨遭遇情況并適當考慮“或”聯(lián)合分布的分析結(jié)果，或基于當前規(guī)范適當增加壩頂高程的安全余度，盡量滿足二者遭遇時的高程要求。