蔡奇霖， 涂金良， 王兆禮， 賴成光

(1.華南理工大學 土木與交通學院， 廣東 廣州 510640； 2.廣東省水利水電技術(shù)中心， 廣東 廣州 510635)

1 研究背景

水庫在防洪、灌溉、供水、發(fā)電等方面對我國國民經(jīng)濟發(fā)展起著不可替代的支撐作用，水庫大壩的安全問題也是涉及人民生命財產(chǎn)安全的一個重大問題[1-2]。據(jù)統(tǒng)計[3]，1954-2006年之間全國有3 498座水庫潰壩。潰壩事故主要發(fā)生在1990年以前，其中1954-1990年共有3 260座水庫潰壩，占總潰壩數(shù)的93.2%，年均約88座；1991-2000年，共有227座水庫潰壩，年均約23 座；2001-2006年，共有35座水庫潰壩，年均約6座。在此期間有2個潰壩高峰，一個是 1960年前后，即1959-1961年間，共計潰壩507座，每年約253座；另一個高峰期在1973年前后，僅1973年就潰壩554座。進入 20 世紀 90年代以來，特別是2001年以后，由于設(shè)計、建設(shè)以及管理水平的不斷提高使得年潰壩數(shù)量明顯減少。雖然如此，每年仍有潰壩事故發(fā)生。毫無疑問，合理的設(shè)計是保證大壩安全運行的前提條件，對保障大壩下游人民的生命財產(chǎn)安全意義重大。

對于水面面積大、周邊空曠的水庫，因其吹程較大容易形成大風浪，導致風壅水位過高以及波浪爬高過大而引起漫頂或躍頂現(xiàn)象，進而威脅大壩的安全。因此，在大壩設(shè)計過程中，確定合理的設(shè)計風速對雍水及風浪爬高計算至關(guān)重要，也決定了大壩的綜合超高以及安全裕度。近幾十年來，全球氣候變暖導致眾多地區(qū)風速和降雨特征均發(fā)生顯著變化[4-6]。不少研究發(fā)現(xiàn)，我國沿海地區(qū)的強風和暴雨事件發(fā)生頻率增加[7-9]，這也意味著兩種事件遭遇的可能性也會進一步提高。與僅考慮單一變量(如降雨或風速)的分析相比，雙變量(如降雨和風速)遭遇的聯(lián)合分布情況更為復雜，這將會為水庫設(shè)計風速的計算帶來極大的挑戰(zhàn)。

對于風速和暴雨聯(lián)合分布方面的研究，國內(nèi)外學者做了大量工作。武占科等[10]對上海地區(qū)的臺風風速和降雨量建立Copula聯(lián)合分布模型，提出了臺風條件下風速和雨強聯(lián)合概率分布函數(shù)，結(jié)果表明Gumbel Copula函數(shù)聯(lián)合概率分布模型適合描述臺風條件下風速和雨強聯(lián)合概率分布。陳立華等[11]對欽州市臺風影響下的風雨聯(lián)合概率分布模型進行了研究，得到了適用于欽州市的風雨聯(lián)合概率分布模型，并對兩場臺風降雨數(shù)據(jù)進行了檢驗，結(jié)果表明該風雨模型結(jié)果可靠，可應用于臺風影響下欽州市降雨情勢的預測分析。Dong等[12]應用Copula函數(shù)對崇明島進行風速和雨強的聯(lián)合回歸概率分析，建立了G-H Copula聯(lián)合概率模型，通過實測數(shù)據(jù)驗證了該模型的效果十分理想。王修勇等[13]以臺風最大風速和最大1 h時段雨量為樣本，選出最優(yōu)Copula風雨聯(lián)合分布函數(shù)模型，并對瓊州海峽大橋風雨極值重現(xiàn)期進行了分析，結(jié)果表明邊緣分布重現(xiàn)期計算出的風速、雨量值均小于聯(lián)合分布重現(xiàn)期的計算值。侯靜惟等[14]以298場熱帶氣旋最大3 s極值風速與總降水量為研究樣本，建立了適用于海南島的熱帶氣旋風雨聯(lián)合概率模型并分析了風雨單個致災因子超過閾值和兩個致災因子超過閾值兩種聯(lián)合重現(xiàn)期，結(jié)果表明聯(lián)合重現(xiàn)期RPand相比于聯(lián)合重現(xiàn)期RPor與歷史熱帶氣旋歸一化損失率具有更好的相關(guān)性。Um等[15]采用Copula函數(shù)分析了風速與臺風降雨的相關(guān)關(guān)系，并將其應用于韓國濟州島的臺風降雨的預測。Bushra等[16]應用Gumbel Copula函數(shù)分析了孟加拉灣沿岸風速與風暴潮的聯(lián)合風險，結(jié)果表明該模型結(jié)果可靠。上述研究均表明，考慮雙變量遭遇情況的計算結(jié)果比僅考慮單一變量更加接近實際情況；然而，上述研究大部分聚焦于相對較大空間尺度的風速和暴雨遭遇問題，而對于沿海地區(qū)水庫在風雨遭遇條件下有關(guān)設(shè)計風速的研究尚不多見。

廣東省沿海地區(qū)屬于東亞季風區(qū)，降水充沛，臺風登陸頻繁。廣東省沿海地區(qū)水庫眾多，但這些水庫無論是在20世紀50、60年代的建設(shè)期間還是在2000年前后的除險加固過程中，其設(shè)計風速均沒有考慮強風暴雨的遭遇影響。在極端暴雨頻次與臺風強度不斷增加的背景下，這些水庫大壩的高程是否仍能滿足安全要求，需要進一步復核。目前的技術(shù)規(guī)范如《碾壓式土石壩設(shè)計規(guī)范》(SL 274—2020)和《砌石壩設(shè)計規(guī)范》(SL 25—2006)是按大壩的重要性給出相應的風速標準，也尚未考慮風速與暴雨的遭遇情況；盡管《小型水利水電工程碾壓式土石壩設(shè)計規(guī)范》(SL 189—2013)要求在沿海地區(qū)應該考慮最大風力與暴雨同時出現(xiàn)的條件，但如何執(zhí)行也沒有明確規(guī)定。

鑒于此，本文將以位于廣東省沿海地區(qū)的東湖水庫為研究對象，探討變化環(huán)境下強風暴雨遭遇的聯(lián)合分布特征，應用Copula函數(shù)構(gòu)建強風-暴雨聯(lián)合概率分布模型，計算設(shè)計風速與暴雨的重現(xiàn)水平，最后基于聯(lián)合分布的設(shè)計風速對水庫主壩壩頂和防浪墻高程進行安全復核。研究成果可為我國沿海地區(qū)水庫大壩的設(shè)計風速計算和壩頂高程設(shè)計提供參考依據(jù)。

2 數(shù)據(jù)來源與研究方法

2.1 數(shù)據(jù)來源

東湖水庫位于廣東省陽江市陽東區(qū)那龍河上游，于1960年建成，2006年經(jīng)加固改造后集雨面積達51.15 km2，水庫總庫容為1.27×108m3，為大(Ⅱ)型水庫，設(shè)計重現(xiàn)期為100年一遇。水庫主壩壩型為均質(zhì)土壩，壩頂高程為36.7 m，壩頂長度為130 m；水庫正常蓄水位為31.00 m，設(shè)計洪水位為32.95 m，校核洪水位為33.73 m。

所采用的降雨和風速數(shù)據(jù)來源于水庫附近的陽江氣象站(東經(jīng)111°58′、北緯21°51′)，其中降雨(逐日)數(shù)據(jù)序列長度為1960-2018年，最大風速(逐日)數(shù)據(jù)序列長度為1971-2018年。提取每年最大日降雨和最大風速作為年最大日降雨和年最大風速，并重新組成年最大日降雨和年最大風速序列；使用滑動法分別獲得年最大3 d和7 d降雨序列，并分別提取二者對應時間段內(nèi)的日最大風速作為年最大風速序列。上述序列均經(jīng)過平穩(wěn)性分析。

2.2 邊緣分布函數(shù)

目前常用的邊緣分布函數(shù)包括皮爾遜Ⅲ型(P-Ⅲ)、伽馬分布(Gamma)、廣義極值分布(generalized extreme value distribution, GEV)、極值Ⅲ型(Weibull)、極值Ⅰ型(Gumbel)等，均取得理想的效果[17]。通過綜合對比和考慮，本研究最終選擇P-Ⅲ、Gamma、GEV 3種常用函數(shù)，各函數(shù)表達式如表1 所示。

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

(1)

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

2.3 Copula函數(shù)

Copula 函數(shù)可描述變量間的相關(guān)性，常用的Copula函數(shù)包括經(jīng)驗Copula函數(shù)、橢圓Copula函數(shù)以及Archimedean Copula函數(shù)等。Archimedean Copula函數(shù)因其計算簡便、形式多樣，在實際工程中得到了廣泛應用[20-23]。本研究將選用Gumbel Copula、Frank Copula 和Clayton Copula 3種常見二元函數(shù)，各Copula函數(shù)表達式見表2。應用AIC(Akaike information criterion)信息準則法[24-25]及OLS(ordinary least squares)離差平方和最小準則[26]進行擬合優(yōu)度評價，AIC值與OLS值越小，則擬合效果越好。

注：表達式中的α、β、γ分別為形狀、尺度、位置參數(shù)。

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

OLS=

(2)

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

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

(3)

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

Copula函數(shù)擬合優(yōu)度與OLS、AIC值呈負相關(guān)關(guān)系。

2.4 聯(lián)合與同現(xiàn)重現(xiàn)期

對水文研究領(lǐng)域的遭遇問題，聯(lián)合風險率和同現(xiàn)風險率是關(guān)注的重點。聯(lián)合風險率表示至少有1個出現(xiàn)非期望事件的概率，同現(xiàn)風險率表示兩個非期望事件同時出現(xiàn)的概率。同現(xiàn)聯(lián)合概率可用公式(4)表示，聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期分別用公式(5)、(6)表示。

P[X1>x1,X2>x2]=1-Fx1(x1)-Fx2(x2)+

F(x1,x2)

(4)

(5)

(6)

式中：P為同現(xiàn)、聯(lián)合概率；X1和X2為變量序列；Fx1和Fx2為邊緣分布函數(shù)；F為聯(lián)合分布函數(shù)；Tor和Tand分別為同現(xiàn)、聯(lián)合重現(xiàn)期。

2.5 壩頂超高計算

選用《碾壓式土石壩設(shè)計規(guī)范》(SL 274—2020)[27]作為壩頂超高計算依據(jù)，其規(guī)定壩頂高程等于水庫靜水位加壩頂超高，應分別按4種運用條件計算并取其最大值，壩頂超高由下式確定：

y=R+e+A

(7)

式中：y為壩頂超高，m；R為波浪在壩坡上的爬高，m；e為風雍水面高度，m；A為安全超高，m。

設(shè)計波浪爬高值應根據(jù)工程等級確定，1、2級壩采用累計頻率為1%的波浪爬高值R1%。由該規(guī)范附錄表A.1.13查得R1%/Rm=2.23。此外，規(guī)范的條件1中建議“采用多年平均最大風速的1.5～2.0倍”，本次分別采用上限2.0倍和下限1.5倍計算，以探討二者計算結(jié)果的差異。

以陽江站1971-2018年間的年最大日、最大3 d、最大7 d 3組降雨序列及其對應年最大風速序列作為計算樣本，應用3種分布函數(shù)和K-S檢驗進行邊緣分布擬合，隨后采用Copula函數(shù)、AIC、OLS準則進行聯(lián)合概率分布擬合，最后計算聯(lián)合與同現(xiàn)重現(xiàn)期下的風速和降雨量。分別采用遭遇風速與規(guī)范推薦的設(shè)計風速計算壩頂高程，并與實際壩頂高程進行對比，綜合復核東湖水庫主壩壩頂高程。

3 結(jié)果與分析

3.1 年最大日降雨量與最大風速變化趨勢

分別對3組降雨及風速序列進行趨勢分析，以年最大日降雨量與年最大風速序列為典型進行分析，二者變化趨勢如圖1、2所示。由圖1、2可知，1960-2018年間，年最大日降雨量呈顯著增大趨勢(P<0.05)；年最大風速總體呈現(xiàn)不顯著增大趨勢(P<0.05)，其中1971-2000年呈減小趨勢而 2001-2018年呈明顯增大趨勢。

圖1 1960-2018年研究區(qū)域年最大日降雨量變化趨勢

3.2 相關(guān)性與遭遇統(tǒng)計

分別計算年最大日、最大3 d與最大7 d降雨量及其對應年最大風速序列的相關(guān)性，可得其相關(guān)系數(shù)分別為0.007、0.128和0.156，表明3組序列相關(guān)性極弱。分別選取年最大風速以及最大日降雨量、最大3 d和最大7 d降雨量統(tǒng)計東湖水庫的風雨遭遇情況(表3)，當東湖水庫發(fā)生年最大風速時剛好遭遇最大日降雨的次數(shù)為2次，占比為4.17%；當發(fā)生年最大風速時恰好遭遇最大3 d降雨的次數(shù)為5次，占比為10.42%；當發(fā)生年最大風速時剛好遭遇最大7 d降雨的次數(shù)為10次，占比為20.83%。以上表明盡管最大風速與最大降雨序列相關(guān)性極弱，但出現(xiàn)了年最大風速和年最大暴雨遭遇的情況。

表3 東湖水庫年最大風速與不同時段最大降雨量遭遇情況統(tǒng)計

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

以1971-2018年間年最大日、最大3 d、最大7 d降雨量及其對應年最大風速為樣本，采用P-Ⅲ、GEV和Gamma 3種平穩(wěn)分布函數(shù)進行邊緣分布擬合，應用K-S、PDE檢驗進行擬合優(yōu)度評價，各序列的邊緣分布擬合優(yōu)度如表4所示。

表4中的擬合優(yōu)度結(jié)果表明，最大日降雨和對應日最大風速的最優(yōu)函數(shù)均為P-Ⅲ分布，最大3 d降雨和對應日最大風速的最優(yōu)函數(shù)分別為GEV分布和P-Ⅲ分布，最大7 d降雨的最優(yōu)函數(shù)為P-Ⅲ分布但對應日最大風速為Gamma分布。

表4 各序列邊緣分布擬合優(yōu)度

采用Gumbel Copula、Frank Copula及Clayton Copula 3種函數(shù)分別構(gòu)建不同序列間的聯(lián)合分布模型，采用AIC信息準則法和OLS準則進行擬合優(yōu)度評價，結(jié)果如表5所示。

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

注：加粗字體表示該邊緣分布函數(shù)擬合優(yōu)度最佳。

注：加粗字體表示該Copula函數(shù)擬合優(yōu)度最佳。

由表5可看出，日最大風速與最大日降雨量聯(lián)合分布的最優(yōu)函數(shù)均為Frank Copula函數(shù)，日最大風速與最大3 d降雨量聯(lián)合分布的最優(yōu)函數(shù)均為Clayton Copula函數(shù)，日最大風速與最大7 d降雨量聯(lián)合分布的最優(yōu)函數(shù)均為Frank Copula函數(shù)，圖3為日最大風速與最大日降雨量、最大3 d降雨量、最大7 d降雨量的聯(lián)合概率分布圖。

圖3 不同序列間聯(lián)合概率分布圖表6 不同重現(xiàn)期和計算工況下壩頂及防浪墻高程計算結(jié)果

3.4 風速與暴雨重現(xiàn)水平

經(jīng)過多種組合計算，發(fā)現(xiàn)強風和暴雨在低重現(xiàn)期遭遇組合條件下對水庫大壩安全性影響有限，因此僅選取二者均為100年一遇(等概率法)遭遇的極端情況的計算結(jié)果進行分析。同現(xiàn)重現(xiàn)期、聯(lián)合重現(xiàn)期以及單變量重現(xiàn)期不同時段下，最大風速和最大降雨量相應壩頂及防浪墻高程計算結(jié)果如表6所示。

由表6可知，基于不同聯(lián)合序列所得出的風速以及降雨量均存在較大差異，其中聯(lián)合重現(xiàn)期計算值最大，同現(xiàn)重現(xiàn)期計算值最小，而單變量重現(xiàn)期位于二者之間。對于最大風速，聯(lián)合重現(xiàn)期比同現(xiàn)重現(xiàn)期大1.5～1.9倍，但聯(lián)合重現(xiàn)期最大風速僅比考慮單變量的最大風速大10%左右；而對于最大降雨量，聯(lián)合重現(xiàn)期則比同現(xiàn)重現(xiàn)期大1.7～1.9倍。

取樣條件同現(xiàn)重現(xiàn)期最大風速/(m·s-1)最大降雨量/mm壩頂高程/m聯(lián)合重現(xiàn)期最大風速/(m·s-1)最大降雨量/mm壩頂高程/m單變量重現(xiàn)期最大風速/(m·s-1)最大降雨量/mm壩頂高程/m2.0倍多年平均最大風速最大風速/(m·s-1)壩頂高程/m1.5倍多年平均最大風速最大風速/(m·s-1)壩頂高程/m實際壩頂高程/m實際防浪墻頂高程/m最大日降雨量14.97393.3435.3823.68681.0836.3021.80616.8636.10最大3 d降雨量16.69557.1135.5627.151106.5036.6724.83960.4736.4235.6737.5926.7536.6336.7037.50最大7 d降雨量16.94653.3635.5833.011221.6237.3129.721106.7936.95

3.5 壩頂高程復核

根據(jù)陽江氣象站實測1971-2018年最大風速資料可求得多年平均最大風速為17.83 m/s。對于遭遇情景，選取同現(xiàn)重現(xiàn)期、聯(lián)合重現(xiàn)期以及單變量重現(xiàn)期的風速最大值并根據(jù)規(guī)范推求壩頂高程。在3種取樣條件下不同重現(xiàn)期壩頂高程最大值、規(guī)范計算值以及實際壩頂高程如表6所示。通過調(diào)查可知，東湖水庫實際壩頂高程為36.70 m，防浪墻頂高程為37.50 m，考慮到東湖水庫大壩防浪墻防御能力較強，因此以防浪墻的高程作為依據(jù)評判是否安全。通過計算發(fā)現(xiàn)，設(shè)計洪水位加正常運行條件下的壩頂超高值最大，因此選取該條件作為壩頂高程計算條件，且分別取規(guī)范的上限2.0倍和下限1.5倍多年平均最大風速進行分析。

結(jié)果表明，當設(shè)計風速取1.5倍多年平均最大風速時，所計算的壩頂高程為36.63 m，當前壩頂和防浪墻實際高程均能滿足要求；而當設(shè)計風速取2.0倍多年平均最大風速時，壩頂高程為37.59 m，比防浪墻實際高程略高。對于遭遇中最危險工況，即聯(lián)合重現(xiàn)期(最大7 d降雨取樣)為100年一遇時所計算的壩頂高程為37.31 m，當前防浪墻高程仍能滿足要求。

4 討 論

不少研究表明[7-9]，最近幾十年來臺風的強度在不斷提高，風速的極值也在不斷增大，而當風速序列增加這些極值樣本時會大大提高多年平均最大風速值。由于東湖水庫在2006年進行了除險加固，當年所使用的風速序列相對較短；由最大風速變化趨勢(圖2)可知在2000年以后年最大風速呈現(xiàn)顯著上升趨勢，但后續(xù)這些極值樣本大部分沒有被納入除險加固的設(shè)計風速計算，這可能是導致當前防浪墻高程略小于2.0倍多年平均最大風速時所推算高程的主要原因。而對于遭遇情況，當考慮強風和暴雨聯(lián)合分布時，同現(xiàn)重現(xiàn)期和聯(lián)合重現(xiàn)期的最大風速與單變量最大風速存在較大差異，其中聯(lián)合重現(xiàn)期最大風速比僅考慮單變量的最大風速大10%左右。通過對比分析可知，東湖水庫的2.0倍多年平均最大風速仍大于聯(lián)合重現(xiàn)期最大風速，表明《碾壓式土石壩設(shè)計規(guī)范》(SL 274—2020)在考慮強風暴雨遭遇時仍可滿足安全要求。此外，通過調(diào)查可知，盡管主壩的防浪墻高程未能滿足2.0倍多年平均最大風速的高程要求，但最近幾十年在臺風登陸前后基本上沒有出現(xiàn)水位過高而導致險情的情況，這主要是因為水庫從未遭遇到2.0倍多年平均最大風速的情景；此外，日益精準的天氣預報以及水庫的科學調(diào)度也有利于水庫提前做好防御措施，例如在獲知臺風預報后可通過迅速泄洪以降低水位，使波高、風浪雍高大大降低，從而保證水庫安全度過險情。當然，本研究僅能說明東湖水庫在考慮強風暴雨遭遇時能滿足安全要求，其它水庫是否滿足仍需要具體分析。盡管通過科學的管理和預報手段可以降低風險，但在全球變暖驅(qū)動下，極端暴雨與臺風強度有可能進一步加強，強風和暴雨聯(lián)合出現(xiàn)的概率也可能會進一步增大，因此在水庫大壩特別是沿海地區(qū)的水庫大壩設(shè)計和建設(shè)過程中，建議增加考慮強風和暴雨遭遇的情況，根據(jù)計算結(jié)果可適當增加壩頂高程的安全裕度，以應對未來發(fā)生極端暴雨和臺風的情況。

圖2 1971-2018年研究區(qū)域年最大風速變化趨勢

5 結(jié) 論

以東湖水庫為研究對象，利用廣義極值分布(GEV)、伽馬分布(Gamma)和皮爾遜Ⅲ型分布(P-Ⅲ)分別對最大風速和降雨序列進行邊緣分布擬合，基于Copula函數(shù)構(gòu)造最大風速和降雨序列的聯(lián)合概率分布，估算100年一遇重現(xiàn)期下最大風速和降雨量設(shè)計值，并復核了水庫主壩壩頂高程，得出如下結(jié)論：

東湖水庫1960-2018年間年最大日降雨量呈顯著上升趨勢(P<0.05)，最大風速則總體呈現(xiàn)不顯著上升趨勢(P<0.05)，其中1971-2000年呈下降趨勢，2001-2018年呈明顯上升趨勢，且出現(xiàn)年最大風速與年最大暴雨遭遇的情況。遭遇分析結(jié)果表明，聯(lián)合重現(xiàn)期計算值最大，同現(xiàn)重現(xiàn)期計算值最小，而單變量重現(xiàn)期位于二者之間。對于最大風速，聯(lián)合重現(xiàn)期比同現(xiàn)重現(xiàn)期大1.5～1.9倍，但聯(lián)合重現(xiàn)期最大風速僅比考慮單變量的最大風速大10%左右。當設(shè)計風速取1.5倍多年平均最大風速時，當前壩頂和防浪墻高程均能滿足要求；當設(shè)計風速取2.0倍多年平均最大風速時，所推算的高程比當前實際防浪墻高程略高。對于遭遇最危險工況(聯(lián)合重現(xiàn)期)，當前的防浪墻高程能滿足要求。東湖水庫的2.0倍多年平均最大風速所計算的防浪墻高程仍大于聯(lián)合重現(xiàn)期最大風速所推算的高程，表明《碾壓式土石壩設(shè)計規(guī)范》(SL 274—2020)在考慮強風暴雨遭遇時仍可滿足要求。但考慮到未來氣候變化，建議增加考慮強風和暴雨遭遇的情況，可適當增加壩頂高程的安全裕度以應對未來發(fā)生極端暴雨和臺風的情況。