王磊之，崔婷婷，李笑天，李伶杰，劉克強，李曦亭

(1.南京水利科學研究院水文水資源與水利工程科學國家重點實驗室，江蘇 南京 210029；2.水利部水資源管理中心，北京 100038； 3.太湖流域管理局水利發(fā)展研究中心，上海 200434)

沿海地區(qū)受海洋性季風的影響，降水較多，加上外海高潮位的頂托作用，強降水、高潮位“碰頭”的情況時有發(fā)生，因此沿海城市的防洪排澇工程系統(tǒng)設(shè)計不僅受區(qū)間暴雨的影響，亦要考慮外海潮位的頂托效應(yīng)。暴雨-潮位的遭遇在重現(xiàn)期意義上直接影響了沿海地區(qū)防洪排澇工程設(shè)計標準，而兩者的遭遇分析通常以其遭遇的聯(lián)合分布為基礎(chǔ)[1]。長久以來，水文學者對兩者遭遇的聯(lián)合重現(xiàn)期開展了不少研究，早期黃國如等[2]基于頻率組合法，以面平均雨量和潮位為研究對象建立了2個因素的水位函數(shù)，推求了感潮河段的設(shè)計洪水位，但其推求的假設(shè)條件為降雨和潮位相互獨立；林榮等[3]采用正態(tài)變換法研究了黃浦江風暴潮位與強降雨遭遇的概率，但該方法在原始數(shù)據(jù)處理過程中會造成部分信息失真。近年來，Copula函數(shù)因其計算過程可逆、推求結(jié)果相對可靠、靈活性和適應(yīng)性較強等優(yōu)勢，在極端水文事件研究，如極端暴雨[4]、洪水分析[5-6]、水文干旱[7]等方面得到廣泛應(yīng)用。在暴雨-潮位遭遇研究方面也出現(xiàn)了不少應(yīng)用Copula函數(shù)計算暴雨-潮位遭遇風險概率的案例。例如：任錦亮等[8]運用Copula函數(shù)構(gòu)建了上海市臺風降水和潮位的聯(lián)合分布模型，對最大1 d降水量與相應(yīng)潮位的遭遇風險概率進行了研究；劉曾美等[9]應(yīng)用Copula函數(shù)計算了中山市中珠聯(lián)圍設(shè)計暴雨和承泄區(qū)潮位的遭遇概率，結(jié)果表明，每年6—7月中珠聯(lián)圍的暴雨有較大概率遭遇西江高潮位；萬永靜等[10]、張衛(wèi)國等[11]分別應(yīng)用Copula函數(shù)分析了南京市、甬江流域設(shè)計極值暴雨與長江最高潮位的遭遇概率；Fu等[12]采用Copula函數(shù)生成降水過程和評估排水系統(tǒng)的效果；劉曾美等[1]、Zheng等[13]、Xu等[14]、徐宗學等[15]分別利用Copula函數(shù)探索了廣東陽山、澳大利亞、?？谑?、深圳市等地沿海地區(qū)暴雨和風暴潮對洪澇災(zāi)害的復(fù)合驅(qū)動作用。

當前對沿海地區(qū)暴雨-潮位遭遇的研究大多集中在采用Copula函數(shù)探討分析長系列極端水文事件要素之間的聯(lián)合分布特征，而對于其聯(lián)合分布特征變異情況的研究卻較為少見。然而，受全球氣候變化的影響，極端水文事件有明顯增加趨勢[16]，吳燕娟[17]通過分析我國不同地區(qū)站點的最大1 d降水量資料，發(fā)現(xiàn)我國東南地區(qū)大雨和暴雨均呈顯著上升趨勢。政府間氣候變化專門委員會(IPCC)第5次評估報告指出，在未來全球氣候變暖的背景下，極端氣象水文事件還將進一步增加[16-18]，因此直接采用長系列數(shù)據(jù)研究暴雨-潮位組合問題無形中忽略了降水特性變化對遭遇特征的影響。本文選擇華南沿海地區(qū)的深圳市為研究區(qū)域，采用Pettitt方法[19]分析檢測極值暴雨的變異情況，對比分析了不同時期深圳市暴雨-潮位遭遇聯(lián)合風險的變化情況，以期為沿海地區(qū)的防洪排澇工作提供參考。

1 研究區(qū)概況與數(shù)據(jù)來源

深圳市地處華南地區(qū)、珠江口東岸，東臨大亞灣和大鵬灣，西瀕珠江口和伶仃洋，南隔深圳河與香港相連。作為中國的經(jīng)濟中心城市和國際化城市，深圳擁有中國最多的口岸數(shù)量和最多的出入境人員，發(fā)達的經(jīng)濟社會以及人口和財富的高度聚集凸顯了深圳市洪澇災(zāi)害的暴露性和脆弱性，加之深圳市瀕臨南海，受季風氣候影響，洪澇災(zāi)害多發(fā)，并易受海洋風暴潮影響，在汛期容易出現(xiàn)極端暴雨、高潮位“碰頭”的情況，使得洪澇災(zāi)害更為嚴重。因此本文選擇深圳市作為典型沿海地區(qū)，重點研究區(qū)域暴雨和外海高潮位遭遇聯(lián)合風險的變化情況。

根據(jù)1965—2016年的《中華人民共和國水文年鑒》，深圳市沿海地區(qū)分布著若干潮位站，其中位于伶仃洋的赤灣站擁有較長系列的潮位及降水數(shù)據(jù)記錄，故本文選擇赤灣站作為研究深圳極值暴雨-高潮位遭遇特征的代表性站點。對深圳市而言，城市防洪除澇對短歷時的最大1 d降水量響應(yīng)程度較高，因此利用赤灣站1965—2016年的逐日降水量、高潮位資料，提取了該站年最大1 d降水量及其對應(yīng)的最高潮位系列。值得注意的是，由于本文研究的是極值降水、高潮位“碰頭”的情況，因此，相應(yīng)的高潮位序列在年最大1 d降水量發(fā)生日期內(nèi)進行提取。

2 研究方法

2.1 突變性診斷

以極值降水(年最大1 d降水量)作為分析降水變化響應(yīng)的特征指標。采用Pettitt方法對深圳赤灣站長系列最大1 d降水量突變性進行診斷，該方法與Mann-Kendall法相似，是一種基于Mann-Whitney統(tǒng)計量的非參數(shù)檢驗方法，其診斷步驟如下：

步驟1構(gòu)造Mann-Kendall秩序列：

(1)

式中：τ為Kendall秩相關(guān)系數(shù)；n為時間序列的長度；x為極值降水序列。

步驟2把所構(gòu)造的秩序列取得最大值S的時刻計為T，認為T為突變時刻，T和S的表達式分別為

T=argmax|Sτ|

(2)

(3)

步驟3計算顯著性統(tǒng)計量P：

P≈2exp[-6S2(n3+n2)]

(4)

若P≤0.5，則認為檢測出來的突變點在統(tǒng)計意義上是顯著的。

2.2 分布擬合與優(yōu)選

2.2.1分布擬合

a.邊緣分布。采用在水文學中常用的5類分布(Normal分布、Lognormal分布、Gamma分布、Logistic分布和Weibull分布)作為赤灣站極值降水和相應(yīng)高潮位系列邊緣分布的備選擬合線型，邊緣分布擬合模型的采納與否通過K-S檢驗的統(tǒng)計量D確定。

b.聯(lián)合分布。Copula函數(shù)是邊緣分布為[0,1]區(qū)間均勻分布的聯(lián)合分布函數(shù)，利用Copula函數(shù)可以將多個隨機變量的邊緣分布連接起來，得到它們的聯(lián)合分布。采用水文學中常用的Archimedean Copula函數(shù)族作為聯(lián)合分布的備用函數(shù)，它是一類重要的Copula函數(shù)，求解簡便，適應(yīng)性強，常用的Archimedean Copula函數(shù)主要包含Clayton Copula、Gumbel Copula和Frank Copula 3類。二維Archimedean Copula函數(shù)的聯(lián)合累積分布、參數(shù)θ和Kendall秩相關(guān)系數(shù)τ的關(guān)系表達式見表1，其中，F(xiàn)(p,z)、C(u,v)均為各聯(lián)合分布的累積分布，前者對應(yīng)于實際的二維變量，p、z分別為2個二維變量；后者對應(yīng)于相應(yīng)的Copula函數(shù)，u、v分別為二維Copula函數(shù)中的2個變量。

表1 二維Archimedean Copula函數(shù)累積分布及θ和τ的關(guān)系

2.2.2分布優(yōu)選

擬合優(yōu)度評價是選擇聯(lián)合概率分布函數(shù)的重要指標[20]，常用的擬合優(yōu)度評價是對經(jīng)驗概率和理論概率進行比較，具體可采用直觀圖形分析法、赤池信息準則(AIC)和離差平方和最小準則(OLS)等[21]進行擬合優(yōu)度評價。本文采用AIC優(yōu)選出極值暴雨-潮位聯(lián)合分布的最優(yōu)Copula函數(shù)，再采用直觀圖形分析法檢驗優(yōu)選出的分布是否能夠滿足擬合精度的要求。

對于二維實測數(shù)據(jù)集，其經(jīng)驗累積頻率由下式計算：

(5)

式中：zi、qj為二維變量的取值；Z、Q為二維變量對應(yīng)的事件；nij為所有實測值對中同時小于或等于(zi,qj)的頻次；N為聯(lián)合實測值總對數(shù)。

AIC為評估、衡量統(tǒng)計模型“擬合”資料優(yōu)良性的一種標準，AIC值計算公式為

I=2k-2lnL

(6)

式中：I為AIC值；k為參數(shù)的數(shù)量；L為似然函數(shù)。I值越小，說明Copula函數(shù)擬合精度越高。

2.3 聯(lián)合風險概率

降水量Z和潮位Q中有一個變量超過某一量級的概率稱為聯(lián)合風險概率，相應(yīng)的重現(xiàn)期稱為聯(lián)合重現(xiàn)期。Z≥z和Q≥q的聯(lián)合風險概率計算公式為

P(Z≥z,Q≥q)=1-F(z,q)

(7)

式中：P(Z≥z,Q≥q)為聯(lián)合風險概率；F(z,q)為Copula函數(shù)計算的聯(lián)合分布概率。

3 結(jié)果與分析

3.1 降水突變分析

圖1為深圳赤灣站1965—2016年年最大1 d降水量序列及其Pettitt統(tǒng)計量的變化趨勢，可見年最大1 d降水量呈現(xiàn)明顯增加趨勢，同時，對年最大1 d降水量的突變檢測表明：Pettitt統(tǒng)計量在1985年取得最小值，因此選擇1985年作為深圳市年最大1 d降水量變化的突變點。這一突變年份與涂新軍等[22-23]在長江流域、珠江流域降水變異情況的研究結(jié)論是較為一致的。

圖1 赤灣站年最大1 d降水量及相應(yīng)Pettitt統(tǒng)計量年際變化

根據(jù)年最大1 d降水量序列的突變情況，將研究時段分為1965—1985年和1986—2016年2個階段，計算2個階段的年最大1 d降水量均值分別為138.2 mm和170.7 mm，第二階段均值明顯大于第一階段均值。

3.2 分布擬合與優(yōu)選

3.2.1邊緣分布

a.分布線型優(yōu)選。采用5類備選邊緣分布線型擬合了赤灣站1965—2016年年最大1d降水量和高潮位系列，并將擬合的邊緣分布結(jié)果采用K-S檢驗和AIC值進行優(yōu)選，結(jié)果見表2。由表2可知，各線型擬合均通過了95%的K-S檢驗，均可作為備選線型；由AIC值可知Lognormal線型是擬合赤灣站極值暴雨的最優(yōu)線型，因此選擇Lognormal分布作為赤灣站極值暴雨的邊緣分布。采用相同方法得出年最大1 d降水量對應(yīng)高潮位的邊緣分布為Normal分布。

表2 赤灣站年最大1 d降水量系列擬合優(yōu)度比較

b.不同階段分布線型變化。依據(jù)分布線型優(yōu)選結(jié)果，采用Lognormal線型、Normal線型分別擬合不同階段赤灣站年最大1 d降水量及其對應(yīng)伶仃洋高潮位的邊緣分布。為直觀描述1965—1985年和1986—2016年2個階段赤灣站降水量、潮位邊緣分布的變化，圖2(a)(b)分別給出了2個階段中降水量、潮位的邊緣分布。由圖2(a)可知，在2個階段中，赤灣站年最大1 d降水量的邊緣分布發(fā)生了顯著變化，如200 mm的年最大1 d降水量在第一階段和第二階段中的對應(yīng)重現(xiàn)期分別約為10 a和3 a，同樣量級的降水量在第二階段中的重現(xiàn)期比在第一階段中要小得多。即同樣量級的降水量在第二階段發(fā)生的概率明顯更高。由圖2(b)可知，赤灣站年最大1 d降水量對應(yīng)高潮位在2個階段的差別要小于年最大1 d降水量，在潮位較高時(高于3.2 m)，同量級潮位在不同階段的重現(xiàn)期相差較小，但在潮位較低時，在第二階段中的重現(xiàn)期要明顯更小?？梢姴煌A段赤灣站最大1 d降水量及其對應(yīng)潮位的邊緣分布均發(fā)生了較大變化，同量級的降水量在第二階段中的重現(xiàn)期均小于在第一階段中的重現(xiàn)期。這說明，深圳市短歷時極值暴雨在新時期內(nèi)均呈現(xiàn)出增大的趨勢。

圖2 赤灣站年最大1 d降水量和相應(yīng)潮位邊緣分布對比

3.2.2聯(lián)合分布

基于不同階段赤灣站年最大1 d降水量及其對應(yīng)高潮位邊緣分布，構(gòu)建赤灣站年最大1 d降水量和潮位的聯(lián)合分布，聯(lián)合分布組合共2個(1965—1985年和1986—2016年)，備選的Archimedean Copula函數(shù)包括Frank Copula、Gumbel Copula和Clayton Copula3類，則共計構(gòu)建6個聯(lián)合分布。聯(lián)合分布類型的優(yōu)選采用AIC法。赤灣站在不同階段年最大1 d降水量和相應(yīng)高潮位遭遇的Copula函數(shù)參數(shù)計算結(jié)果如表3所示。

表3 赤灣站年最大1 d降水量及相應(yīng)潮位系列

利用表3中的Copula函數(shù)參數(shù)分別構(gòu)建了2個階段赤灣站年最大1 d降水量和對應(yīng)最高潮位的聯(lián)合分布，并采用各Copula函數(shù)的累積概率公式分別計算了歷年最大1 d降水量及對應(yīng)潮位組合的理論概率。圖3為3類Copula函數(shù)計算的赤灣站年最大1 d降水量和對應(yīng)高潮位經(jīng)驗頻率和理論頻率的對比，3類Copula函數(shù)擬合的理論-經(jīng)驗數(shù)據(jù)點基本集中在y=x附近，可見3類Copula函數(shù)對深圳年最大1 d降水量與對應(yīng)高潮位的擬合效果均較好。

圖3 赤灣站年最大1 d降水量與對應(yīng)潮位Copula經(jīng)驗頻率和理論頻率對比Fig.3 Comparison of experimental and theoretical frequencies of annual maximum one-day precipitation and corresponding tidal level at Chiwan Station

采用式(6)計算了各Copula函數(shù)擬合暴雨-潮位聯(lián)合分布的經(jīng)驗-理論頻率的AIC值，F(xiàn)rank Copula、Clayton Copula、Gumbel Copula擬合的AIC值分別為-176.3、-172.3和-168.3，可知Frank Copula擬合的精度略高于其余2類Copula函數(shù)。因此，本文采用基于Frank Copula函數(shù)的聯(lián)合分布模型構(gòu)建赤灣站年最大1 d降水量-高潮位聯(lián)合分布。

將赤灣站暴雨-潮位區(qū)間劃分為若干網(wǎng)格，計算了2個階段聯(lián)合分布的概率密度和累積分布，結(jié)果分別如圖4和圖5所示。由圖4和圖5可知，同樣的暴雨-潮位組合在2個階段中的概率密度和累積分布具有明顯不同，總體來說，暴雨-潮位遭遇在第二階段的概率密度要超過第一階段，這導致其在第二階段的累積分布也超過第一階段。從圖5進一步可以看出，年最大1 d降水量的變化對累積分布曲面施加的影響明顯更加強烈，而相應(yīng)高潮位的變化對累積分布曲面施加的影響相對較小，這一表現(xiàn)與前文年最大1 d降水量以及潮位邊緣分布的變化情況是一致的。

(a) 1965—1985年

(a) 1965—1985年

3.3 暴雨-潮位遭遇聯(lián)合風險分析

年最大1 d降水量選擇50 a、20 a重現(xiàn)期，相應(yīng)的高潮位選擇50 a、20 a、10 a、5 a、2 a重現(xiàn)期，在兩兩組合的情景下分別計算了赤灣站在不同階段年最大1 d降水量與相應(yīng)潮位遭遇的聯(lián)合風險，并給出了第二階段風險值相對于第一階段的變化率，結(jié)果見表4。值得注意的是，同樣重現(xiàn)期的年最大1 d降水量及相應(yīng)潮位是指在全系列(1965—2016年)中的重現(xiàn)期，因此在前后2個階段是統(tǒng)一的。由表4可知：①1986—2016年不同的暴雨-潮位遭遇的聯(lián)合風險高于1965—1985年，如對于5 a一遇年最大1 d降水量，聯(lián)合風險增長率為71%～128%，這意味著同樣規(guī)模的暴雨-潮位組合在第二階段的聯(lián)合風險要遠超第一階段；②高潮位重現(xiàn)期越大的暴雨-潮位組合，第二階段聯(lián)合風險的增長率相對越小，如對于50 a一遇年最大1 d降水量，隨著極值暴雨遭遇的潮位重現(xiàn)期從50 a一遇降低為2 a一遇，第二階段相對于第一階段聯(lián)合風險的增長率相應(yīng)從71%增大至128%。

表4 赤灣站年最大1 d降水量與相應(yīng)潮位遭遇聯(lián)合風險

4 結(jié) 論

a.1985年為深圳市年最大1 d降水序列的突變年份，以該年份為依據(jù)劃分研究時段為前后2個階段(1965—1985年和1986—2016年)。經(jīng)優(yōu)選，確定Lognormal線型、Normal線型分別為不同階段赤灣站年最大1 d降水量、相應(yīng)高潮位的最優(yōu)邊緣分布，而Frank Copula為赤灣站年最大1 d降水量-高潮位聯(lián)合分布的最優(yōu)聯(lián)合分布函數(shù)。

b.當年最大1 d降水量重現(xiàn)期為50 a一遇和20 a一遇時，1985年后年最大1 d降水量和不同重現(xiàn)期相應(yīng)高潮位遭遇的聯(lián)合風險相對于1985年前增加了71%～128%，這意味著同樣的暴雨-潮位組合在第二階段的聯(lián)合風險遠超第一階段，這一風險的顯著增加值得水利和市政部門的關(guān)注。

c.Copula函數(shù)作為一類求解簡便、適應(yīng)性強的聯(lián)合分布分析工具，可用于定量分析深圳市暴雨-潮位遭遇的聯(lián)合風險的變化情況，也可為其他沿海城市和地區(qū)的防洪排澇工作提供參考。