基于貝葉斯推斷的長(zhǎng)江口復(fù)合洪水潛在風(fēng)險(xiǎn)

關(guān)鍵詞：復(fù)合洪水；風(fēng)暴潮；不確定性；貝葉斯推斷；Copula方法；長(zhǎng)江口

河口三角洲是洪澇災(zāi)害最活躍的地區(qū)之一[1]。20世紀(jì)70年代以來，全球海洋持續(xù)變暖，亞洲東部沿海地區(qū)登陸的強(qiáng)臺(tái)風(fēng)數(shù)量和伴隨的降水量明顯增加[2-3]。河口地區(qū)通常地勢(shì)低洼，受徑流和天文潮雙相動(dòng)力控制，在臺(tái)風(fēng)登陸時(shí)還會(huì)受風(fēng)暴潮和強(qiáng)降雨的影響，造成復(fù)合洪水災(zāi)害。復(fù)合洪水是指2個(gè)以上的災(zāi)害因子同時(shí)達(dá)到閾值或相隔時(shí)間較短導(dǎo)致洪水災(zāi)害后果急劇放大的現(xiàn)象[4]。與單一驅(qū)動(dòng)因子的極端事件相比，復(fù)合極端事件的潛在破壞性更大[5]。臺(tái)風(fēng)登陸時(shí)不僅會(huì)抬高沿海水位，還會(huì)伴隨大范圍、高強(qiáng)度的降雨，導(dǎo)致內(nèi)河流量激增水位暴漲。外海潮位的抬升又會(huì)阻塞河道洪水下泄，雍高河道水位[6]，導(dǎo)致沿海城市重力排水系統(tǒng)無法正常工作，造成城市雨洪。在氣候變化背景下，評(píng)估河口復(fù)合洪水的潛在風(fēng)險(xiǎn)，對(duì)于防范和應(yīng)對(duì)河口極端洪澇事件具有重要意義[7]。

近年來，海岸復(fù)合洪水災(zāi)害的頻繁發(fā)生引起政府和科學(xué)界的廣泛關(guān)注[8]。多個(gè)時(shí)空尺度的研究表明，沿海洪水（或風(fēng)暴增水）的頻率或強(qiáng)度與降雨（或徑流）在統(tǒng)計(jì)上存在著相關(guān)性[4，9-10]。河口地區(qū)復(fù)合洪水災(zāi)害一直存在被嚴(yán)重低估的現(xiàn)象[11]。Marcos等[12]發(fā)現(xiàn)在全球55%的海岸線極端風(fēng)暴潮和臺(tái)風(fēng)浪會(huì)同時(shí)出現(xiàn)，如果忽略兩者相關(guān)性，30%海岸的極端海平面重現(xiàn)風(fēng)險(xiǎn)將低估50%以上；Kew等[13]發(fā)現(xiàn)荷蘭萊茵河20日極端累計(jì)降雨后出現(xiàn)極值風(fēng)暴潮的概率是將兩者視為相互獨(dú)立的3倍。在實(shí)際工程中，適當(dāng)提高防洪水位設(shè)計(jì)值可有效降低復(fù)合洪水災(zāi)害頻率[14]。Santos等[15]發(fā)現(xiàn)當(dāng)考慮驅(qū)動(dòng)因子之間（風(fēng)暴潮和徑流）的相關(guān)性時(shí)，美國得克薩斯州原防洪設(shè)計(jì)頻率的1%水位比實(shí)際值低0.35m。臺(tái)風(fēng)哈維（Harvey）期間發(fā)生上游河道洪峰、天文大潮和風(fēng)暴潮遭遇的“三碰頭”現(xiàn)象，致使美國得克薩斯州東南部沿海實(shí)際風(fēng)暴潮比預(yù)測(cè)值高3m以上[16]。Sebastian等[17]、Couasnon等[18]采用貝葉斯方法生成復(fù)合洪水事件并進(jìn)行了大量數(shù)值模擬，研究結(jié)果強(qiáng)調(diào)了災(zāi)害因子依賴性對(duì)洪水風(fēng)險(xiǎn)評(píng)估的重要性；Naseri等[19]采用基于貝葉斯的非平穩(wěn)框架對(duì)復(fù)合洪水進(jìn)行了研究，發(fā)現(xiàn)與海平面和依賴性相比，降水對(duì)不確定性的貢獻(xiàn)最大。可見，隨著氣候變化的加劇，僅考慮單一變量超越風(fēng)險(xiǎn)的設(shè)計(jì)標(biāo)準(zhǔn)逐漸難以滿足沿海地區(qū)的防洪安全需求。在中國沿海地區(qū)，70%以上造成災(zāi)害的風(fēng)暴潮均為復(fù)合洪水[20]。長(zhǎng)江口是中國第一大河口，風(fēng)暴潮“二碰頭”“三碰頭”屢屢發(fā)生[21-22]。然而，現(xiàn)有文獻(xiàn)多基于物理過程（如數(shù)值模擬等）研究臺(tái)風(fēng)期間長(zhǎng)江口潮位變化[23]，尚缺乏對(duì)長(zhǎng)江口地區(qū)復(fù)合洪水風(fēng)險(xiǎn)的具體量化和不確定性分析等相關(guān)研究。

本文基于長(zhǎng)周期實(shí)測(cè)數(shù)據(jù)分析長(zhǎng)江口極值水位變化趨勢(shì)及其成因，結(jié)合近岸高分辨率海洋數(shù)值模型重構(gòu)長(zhǎng)江口復(fù)合洪水災(zāi)害因子數(shù)據(jù)集，采用基于貝葉斯框架的Copula函數(shù)對(duì)風(fēng)暴潮“二碰頭”“三碰頭”災(zāi)害潛在風(fēng)險(xiǎn)進(jìn)行研究，并對(duì)復(fù)合洪水風(fēng)險(xiǎn)的不確定性進(jìn)行分解。

1資料和方法

1.1數(shù)據(jù)說明

（1）數(shù)據(jù)介紹。為研究長(zhǎng)江口復(fù)合洪水災(zāi)害，本文收集了1959—2020年降雨、徑流、天文潮和風(fēng)暴潮等4類水文資料。其中，降雨資料選用上海寶山氣象站日累計(jì)降水量；徑流資料選用大通水文站逐時(shí)流量資料；長(zhǎng)江口內(nèi)潮位資料僅收集到徐六涇水文站1988—2020年逐時(shí)潮位。為補(bǔ)全1959—1987年的潮位資料空缺，采用FVCOM（Finite-VolumeCommunityOceanModel）構(gòu)建的三維高分辨率數(shù)值模型，對(duì)長(zhǎng)江口早期的天文潮和風(fēng)暴潮進(jìn)行了批量數(shù)值模擬。除此之外，還收集了1988—2020年長(zhǎng)江口外海大戢山潮位站海平面變化數(shù)據(jù)。

（2）數(shù)據(jù)預(yù)處理。本文選取了1959—2020年發(fā)生在長(zhǎng)江口的148場(chǎng)臺(tái)風(fēng)，并對(duì)該期間降雨、徑流、天文高潮位、風(fēng)暴潮高潮位和增水進(jìn)行了計(jì)算和提取。具體方法為：①提取148場(chǎng)臺(tái)風(fēng)期間的最高潮位，即風(fēng)暴潮高潮位；②通過調(diào)和分析和數(shù)值模擬計(jì)算臺(tái)風(fēng)期間對(duì)應(yīng)的天文潮位，并提取風(fēng)暴潮高潮位同一潮周期內(nèi)的天文潮高潮位；③提取風(fēng)暴增水，即風(fēng)暴潮高潮位和天文潮高潮位之差；④提取風(fēng)暴潮高潮位出現(xiàn)當(dāng)日的累計(jì)降水量；⑤提取風(fēng)暴潮高潮位出現(xiàn)前6日的大通流量作為上游洪水?dāng)?shù)據(jù)。通過上述方法重構(gòu)了包括降雨、徑流、天文潮高潮位、風(fēng)暴潮高潮位和增水的長(zhǎng)江口復(fù)合洪水5種災(zāi)害因子數(shù)據(jù)，其維度為148×5。

1.2方法簡(jiǎn)介

1.2.1基于貝葉斯推斷的Copula方法

采用馬爾可夫鏈蒙特卡洛（MCMC）算法對(duì)Copula參數(shù)進(jìn)行估計(jì)。MCMC算法是一種基于貝葉斯框架的采樣方法，其利用馬爾可夫鏈來構(gòu)建一個(gè)隨機(jī)游走過程，最終收斂到參數(shù)的后驗(yàn)分布。在貝葉斯框架中，可以根據(jù)觀測(cè)數(shù)據(jù)和先驗(yàn)分布計(jì)算單變量或雙變量模型參數(shù)的后驗(yàn)分布，尤其適用于觀測(cè)數(shù)據(jù)較少的情況。與傳統(tǒng)最大似然估計(jì)相比，MCMC算法不僅能夠估計(jì)參數(shù)值，還能量化參數(shù)的不確定性。貝葉斯定理可表示為

式中：為后驗(yàn)分布，即待估計(jì)的模型參數(shù)分布；P（φ）為先驗(yàn)分布；P（Yobsjφ）為似然函數(shù)；P（Yobs）為量綱一化常數(shù)。

對(duì)于二元Copula函數(shù)的參數(shù)求解，似然函數(shù)可按下式計(jì)算：

式中：n為數(shù)據(jù)長(zhǎng)度；C為參數(shù)集為 的Copula函數(shù)；F（ ）和G（ ）分別為變量X和Y的邊緣分布。在確定了最佳邊緣分布、Copula函數(shù)類型和先驗(yàn)分布后，對(duì)式（2）取對(duì)數(shù)后代入式（1），則Copula函數(shù)參數(shù)的后驗(yàn)分布可表示為

1.2.2聯(lián)合分布的重現(xiàn)期

重現(xiàn)期指超過某一閾值的事件連續(xù)2次發(fā)生的平均時(shí)間間隔[24]。在二元洪水風(fēng)險(xiǎn)研究中，同現(xiàn)重現(xiàn)期和聯(lián)合重現(xiàn)期（又稱AND重現(xiàn)期和OR重現(xiàn)期）得到廣泛應(yīng)用[25]。AND重現(xiàn)期描述了2個(gè)邊緣變量同時(shí)超過特定閾值事件發(fā)生的風(fēng)險(xiǎn)性，而OR重現(xiàn)期描述了至少1個(gè)邊緣變量超過特定閾值事件發(fā)生的風(fēng)險(xiǎn)性，兩者公式如下：

1.2.3不確定性及其分解

（1）自然變化導(dǎo)致的重現(xiàn)期變化范圍。不確定性研究可以更準(zhǔn)確地評(píng)估極端事件的風(fēng)險(xiǎn)。根據(jù)Bevacqua等[26]的研究，自然變化導(dǎo)致的重現(xiàn)期變化范圍主要通過Bootstrap抽樣進(jìn)行計(jì)算。具體方法如下：以降雨-潮位二元變量數(shù)據(jù)為例，有放回的抽次降雨潮位數(shù)據(jù)對(duì)（為原始數(shù)據(jù)對(duì)長(zhǎng)度），組成新的序列后計(jì)算其重現(xiàn)期。將上述過程重復(fù)次即可得到包含個(gè)重現(xiàn)期的新數(shù)據(jù)集。該數(shù)據(jù)集的2.5th和97.5th百分位數(shù)（）則為95%置信區(qū)間。相應(yīng)的重現(xiàn)期變化的95%置信區(qū)間則為原始數(shù)據(jù)計(jì)算的重現(xiàn)期T與（）之間的百分比差異，即（）/T。

（2）不確定性的分解。對(duì)于二維聯(lián)合分布，復(fù)合洪水風(fēng)險(xiǎn)評(píng)估結(jié)果不確定性的來源有二：邊緣分布的不確定性和相關(guān)性的不確定性。為了量化單個(gè)因子對(duì)不確定性的貢獻(xiàn)，采用MCMC算法采樣獲取降雨、潮位以及Copula參數(shù)后驗(yàn)分布各5000組數(shù)據(jù)。當(dāng)評(píng)估降雨對(duì)不確定性的貢獻(xiàn)時(shí)，則允許降雨邊緣分布參數(shù)的變化，潮位邊緣分布和Copula參數(shù)取各自分布的中位數(shù)。通過式（4）或式（5）可以計(jì)算出5000組重現(xiàn)期，可得重現(xiàn)期的第5、50和95百分位數(shù)（、和），則降雨導(dǎo)致的重現(xiàn)期的不確定性對(duì)稱范圍（）計(jì)算公式如下[19]：

同樣地，可以計(jì)算出潮位和Copula參數(shù)導(dǎo)致的重現(xiàn)期不確定性對(duì)稱范圍。則第個(gè)變量對(duì)重現(xiàn)期不確定性的相對(duì)貢獻(xiàn)（）可按下式計(jì)算：

需要注意的是，重現(xiàn)期并非單個(gè)變量的線性疊加，因而降雨、潮位和Copula函數(shù)三者的ΔTSym之和并不等于總的重現(xiàn)期不確定性[26]。

2結(jié)果及分析

2.1長(zhǎng)江口極值水位變化趨勢(shì)及成因

長(zhǎng)江口年特征水位變化趨勢(shì)如圖1所示?？梢钥闯觯ǔＧ闆r下年最高水位不超過6m（吳淞基面）。歷年最高潮位為6.73m，出現(xiàn)在1997年的臺(tái)風(fēng)溫妮（No.9711）期間。1988—2020年長(zhǎng)江口徐六涇站的年最高水位整體呈現(xiàn)出明顯的下降趨勢(shì)，下降速率為11.6mm/a。徐六涇站附近的洪水強(qiáng)度雖有所下降，但頻繁發(fā)生的極端水位事件（如1997年、2000年和2002年）表明長(zhǎng)江口的洪水風(fēng)險(xiǎn)依然不容忽視。與外海相對(duì)穩(wěn)定的海平面不同，徐六涇站汛期平均水位的年際變化較為劇烈，并且同樣呈現(xiàn)出波動(dòng)下降的趨勢(shì)，下降速率為1.6mm/a，但兩者相關(guān)性較弱，僅為0.292。徐六涇站的汛期平均水位不僅受到全球海平面上升趨勢(shì)的影響，還受上游徑流的顯著影響。從圖1可以看出，1988—2020年長(zhǎng)江大通汛期平均流量呈現(xiàn)出年際下降趨勢(shì)，下降速率為124.4m3/s。相較于外海海平面，徐六涇站的汛期平均水位與同期上游流量的變化關(guān)系更為密切，即當(dāng)汛期大通流量增加時(shí)，徐六涇站的平均水位較上一年有所抬升；而當(dāng)大通流量減少時(shí)，平均水位則相應(yīng)下降，兩者的相關(guān)性（rtw）高達(dá)0.886。在氣候變暖背景下，大部分海域的極值潮位都呈現(xiàn)出上升趨勢(shì)。然而，過濾掉平均海面變化后這種趨勢(shì)顯著減緩，表明極值水位上升主要由海平面上升驅(qū)動(dòng)[27-28]。與此不同的是，徐六涇站年最高潮位的下降是上游徑流減少和風(fēng)暴潮強(qiáng)度減弱共同作用的結(jié)果，其中風(fēng)暴潮變化是主導(dǎo)因素。

2.2復(fù)合洪水邊緣分布的擬合與檢驗(yàn)

選用6種常見的水文變量分布，即正態(tài)分布、指數(shù)分布、Logistic分布、伽馬分布、PearsonⅢ型分布和廣義極值分布（GeneralizedExtremeValueDistribution，GEV），分別對(duì)各變量時(shí)間序列進(jìn)行擬合（各分布表達(dá)式見文獻(xiàn)[29]），并采用K-S檢驗(yàn)、經(jīng)驗(yàn)分布和核密度估計(jì)圖對(duì)擬合結(jié)果進(jìn)行檢驗(yàn)。表1為最佳邊緣分布擬合和K-S檢驗(yàn)統(tǒng)計(jì)表?？梢钥闯?，大通徑流量分布類型為L(zhǎng)ogistic分布；日累計(jì)降水量分布類型均為Gamma分布；增水和風(fēng)暴潮高潮位分布類型均為PearsonⅢ分布；天文潮高潮位分布類型為GEV分布。各災(zāi)害因子邊緣分布函數(shù)的理論分布與經(jīng)驗(yàn)分布的均方根誤差（ERMS）值基本都在0.035以內(nèi)，且均通過了K-S檢驗(yàn)（P＞0.05）。

圖2為各災(zāi)害因子理論分布、經(jīng)驗(yàn)分布與核密度估計(jì)對(duì)比圖?？梢钥闯?，各災(zāi)害因子的理論分布、經(jīng)驗(yàn)分布和核密度估計(jì)曲線符合較好，且核密度估計(jì)曲線較為平滑，計(jì)算較為復(fù)雜時(shí)可用核密度估計(jì)值代替理論分布值。因此，所選邊緣分布較為合理，可用于多元聯(lián)合分布相關(guān)計(jì)算。

2.3風(fēng)暴潮“二碰頭”災(zāi)害重現(xiàn)期分析和不確定性分解

風(fēng)暴潮“二碰頭”復(fù)合洪水災(zāi)害風(fēng)險(xiǎn)性主要表現(xiàn)在兩方面：一是臺(tái)風(fēng)期間風(fēng)暴潮會(huì)抬升沿海潮位，威脅海堤及近海低洼地區(qū)安全；二是臺(tái)風(fēng)登陸前后一般會(huì)伴隨大范圍、高強(qiáng)度的暴雨，沿海水位的抬升會(huì)導(dǎo)致雨洪難以排出入海，增加濱海城市內(nèi)澇的風(fēng)險(xiǎn)。如2022年6月15—16日期間，珠江流域的強(qiáng)降雨僅與天文大潮高潮位遭遇，便造成流域內(nèi)100余條河流發(fā)生超預(yù)警洪水，部分地區(qū)發(fā)生流域性大洪水。在臺(tái)風(fēng)暴雨期間，海面高度的抬升為洪水提供了更高的基礎(chǔ)水位，嚴(yán)重削弱了河口城市的防洪排洪能力。在降雨-風(fēng)暴高潮位兩災(zāi)害因子的研究中選用Gaussian、Studentt、Clayton、Gumbel、Frank和Joe6類常用Copula函數(shù)作為備選的Copula函數(shù)，其估計(jì)值和檢驗(yàn)結(jié)果見表2?？梢钥闯?，6個(gè)Copula函數(shù)均通過了擬合優(yōu)度檢驗(yàn)，但從AIC統(tǒng)計(jì)量結(jié)果看，GumbelCopula更適合描述臺(tái)風(fēng)降雨-風(fēng)暴高潮位之間的聯(lián)合分布關(guān)系。GumbelCopula函數(shù)是極值Copula的一種，主要用以描述正相關(guān)的（上尾）極值關(guān)系。

圖3為基于GumbelCopula函數(shù)的風(fēng)暴高潮位和臺(tái)風(fēng)暴雨聯(lián)合分布的密度函數(shù)與分布函數(shù)。從圖3中可以看出，當(dāng)風(fēng)暴高潮位抬升時(shí)，臺(tái)風(fēng)降雨強(qiáng)度也隨之增大，這表明在長(zhǎng)江口地區(qū)臺(tái)風(fēng)登陸和發(fā)生極端風(fēng)暴潮高潮位的同時(shí)，遭遇極端降雨的可能性非常大。

圖4和圖5分別為臺(tái)風(fēng)暴雨與風(fēng)暴潮高潮位聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期分布圖。臺(tái)風(fēng)暴雨日累計(jì)降水量的10a、20a、50a、100a重現(xiàn)期對(duì)應(yīng)的重現(xiàn)水平分別為92.2、122.5、164.4和197.0mm；風(fēng)暴潮高潮位相應(yīng)重現(xiàn)期的重現(xiàn)水平分別為5.89、6.15、6.46和6.68m。根據(jù)聯(lián)合分布函數(shù)，10年一遇的臺(tái)風(fēng)暴雨遭遇10年一遇的風(fēng)暴潮高潮位的聯(lián)合重現(xiàn)期為6a，同現(xiàn)重現(xiàn)期為42a；20年一遇的兩者聯(lián)合重現(xiàn)期為11a，同現(xiàn)重現(xiàn)期為89a；50年一遇的兩者聯(lián)合重現(xiàn)期為28a，同現(xiàn)重現(xiàn)期為232a；100年一遇的兩者聯(lián)合重現(xiàn)期為56a，同現(xiàn)重現(xiàn)期為471a?？梢钥闯觯?lián)合重現(xiàn)期的計(jì)算值偏小，2個(gè)100年一遇的變量的聯(lián)合重現(xiàn)期約為50a，采用聯(lián)合重現(xiàn)期作為設(shè)計(jì)值是偏保守的；同現(xiàn)重現(xiàn)期的計(jì)算值偏大，2個(gè)100年一遇的變量的同現(xiàn)重現(xiàn)期約為500a，因而采用同現(xiàn)重現(xiàn)期計(jì)算的設(shè)計(jì)值是稍偏激進(jìn)的。而當(dāng)不考慮變量之間的相關(guān)性時(shí)，可能造成復(fù)合洪水災(zāi)害風(fēng)險(xiǎn)的低估。以同現(xiàn)重現(xiàn)期為參考（圖6），當(dāng)忽略降雨和風(fēng)暴潮高潮位之間的相關(guān)性時(shí)，兩者10年一遇重現(xiàn)期的復(fù)合洪水風(fēng)險(xiǎn)將低估58.2%；20年一遇重現(xiàn)期的復(fù)合洪水風(fēng)險(xiǎn)將低估77.7%；50年一遇重現(xiàn)期的復(fù)合洪水風(fēng)險(xiǎn)將低估90.7%；100年一遇重現(xiàn)期的復(fù)合洪水風(fēng)險(xiǎn)將低估95.3%?？梢姡?dāng)不考慮變量之間相關(guān)性時(shí)，災(zāi)害因子的重現(xiàn)期越長(zhǎng)，復(fù)合洪水災(zāi)害重現(xiàn)期的低估程度越高。

長(zhǎng)江口地區(qū)在臺(tái)風(fēng)登陸和發(fā)生極端風(fēng)暴潮高潮位的同時(shí)，遭遇極端降雨的可能性較大。因此，河口城市在提高沿海堤防防護(hù)標(biāo)準(zhǔn)的同時(shí)，也要考慮極端潮位與強(qiáng)降雨遭遇的風(fēng)險(xiǎn)，增強(qiáng)城市雨洪的排泄能力，如合理布設(shè)城市內(nèi)雨洪行進(jìn)通道，增大城市的人工排澇能力，最大程度上減輕風(fēng)暴潮引發(fā)的水位抬升對(duì)城市排水系統(tǒng)的影響，避免自然排水受阻導(dǎo)致的城市內(nèi)澇風(fēng)險(xiǎn)。

復(fù)合災(zāi)害風(fēng)險(xiǎn)具有一定的不確定性。當(dāng)不確定性較大時(shí)，復(fù)合洪水風(fēng)險(xiǎn)評(píng)估的可信度難以保證。不確定性分解是提升風(fēng)險(xiǎn)評(píng)估可信度和優(yōu)化決策分析的重要步驟。通過分解不確定性，可以識(shí)別對(duì)結(jié)果影響最大的因素，從而有針對(duì)性地進(jìn)行改進(jìn)。以長(zhǎng)江口風(fēng)暴潮“二碰頭”復(fù)合洪水事件為例，根據(jù)Bevacqua等[26]的研究，設(shè)定Bootstrap法抽樣次數(shù)為700。計(jì)算結(jié)果顯示，自然變化導(dǎo)致的重現(xiàn)期變化的95%置信區(qū)間為[?12.6%，16.3%]，即自然變化導(dǎo)致的重現(xiàn)期不確定性為[?12.6%，16.3%]。對(duì)于二維聯(lián)合分布，復(fù)合洪水風(fēng)險(xiǎn)評(píng)估結(jié)果不確定性的主要來源是邊緣分布的不確定性和相關(guān)性的不確定性。為了量化各單一因子對(duì)不確定性的貢獻(xiàn)，采用MCMC方法和式（7）分別對(duì)降雨、風(fēng)暴潮高潮位以及Copula參數(shù)進(jìn)行5000次抽樣并計(jì)算其重現(xiàn)期[19]。結(jié)果顯示，降雨、風(fēng)暴潮高潮位和依賴性對(duì)復(fù)合洪水災(zāi)害不確定性的貢獻(xiàn)分別為54%、41%和5%?？梢姡L(zhǎng)江口復(fù)合洪水風(fēng)險(xiǎn)不確定性的主要來源是降雨。

2.4風(fēng)暴潮“三碰頭”災(zāi)害重現(xiàn)期分析

河口地區(qū)受徑流和潮汐雙相動(dòng)力的控制。如果河口風(fēng)暴潮發(fā)生在天文大潮期間，同時(shí)遭遇上游洪水，將會(huì)帶來巨大的淹沒風(fēng)險(xiǎn)，這種現(xiàn)象被稱為風(fēng)暴潮“三碰頭”現(xiàn)象。單一參數(shù)的Copula函數(shù)可以很好地描述兩變量的聯(lián)合分布。但對(duì)于更高維分布（三變量及以上），各變量之間的相互依賴關(guān)系往往呈現(xiàn)出非對(duì)稱性，因此非對(duì)稱的嵌套Copula更加適用。目前，非對(duì)稱的嵌套Copula已在洪水和干旱特征等方面得到一些應(yīng)用[29]，但在風(fēng)暴潮“三碰頭”災(zāi)害風(fēng)險(xiǎn)量化方面的研究相對(duì)較少。本節(jié)采用非對(duì)稱嵌套Copula對(duì)風(fēng)暴潮“三碰頭”進(jìn)行建模，其結(jié)構(gòu)如圖7，將天文高潮位-風(fēng)暴增水關(guān)系設(shè)置為內(nèi)層Copula（GumbelCopula），大通徑流與內(nèi)層Copula結(jié)果設(shè)置為外層Copula（SurvivalJoeCopula，參數(shù)為），且通過擬合優(yōu)度檢驗(yàn)。結(jié)果表明，大通流量和臺(tái)風(fēng)及其伴隨事件的相關(guān)性較弱。圖8為臺(tái)風(fēng)登陸（靠近）長(zhǎng)江口期間大通流量變異系數(shù)圖，圖中黑色圓圈為每場(chǎng)臺(tái)風(fēng)靠近長(zhǎng)江口期間最高風(fēng)暴潮位前7日內(nèi)大通日流量的變異系數(shù)（CV，方差/均值），藍(lán)色為變異系數(shù)的小提琴圖?？梢钥闯觯儺愊禂?shù)的范圍一般在0～0.1之間，其中25%分位數(shù)為0.016，88%分位數(shù)為0.050，中位數(shù)為0.024?？傮w來看，臺(tái)風(fēng)在靠近長(zhǎng)江口期間，大通日流量變化一般不超過10%，且多數(shù)情況下在5%以內(nèi)，可見臺(tái)風(fēng)系統(tǒng)對(duì)大通日流量的影響較小。因此，在臺(tái)風(fēng)及伴生事件與大通流量的復(fù)合洪水事件研究中，亦可將兩者視為獨(dú)立。根據(jù)三維聯(lián)合分布函數(shù)，計(jì)算長(zhǎng)江口徐六涇水文站徑流、天文高潮位和風(fēng)暴增水“三碰頭”復(fù)合洪水災(zāi)害風(fēng)險(xiǎn)率。對(duì)于聯(lián)合風(fēng)險(xiǎn)率：10年一遇的徑流、天文潮和風(fēng)暴潮增水相遇的聯(lián)合重現(xiàn)期為4a，20年一遇的三者聯(lián)合重現(xiàn)期為7a，50年一遇的三者聯(lián)合重現(xiàn)期為18a，100年一遇的三者聯(lián)合重現(xiàn)期為35a。

3結(jié)論

本文以長(zhǎng)江口為研究區(qū)域，結(jié)合實(shí)測(cè)資料和數(shù)值模擬構(gòu)建了1959—2020年復(fù)合洪水災(zāi)害因子數(shù)據(jù)集，采用基于貝葉斯推斷的Copula相關(guān)方法對(duì)長(zhǎng)江口風(fēng)暴潮“二碰頭”“三碰頭”復(fù)合災(zāi)害風(fēng)險(xiǎn)性進(jìn)行了研究，并對(duì)不確定性進(jìn)行了分析和分解。主要結(jié)論如下：

（1）1959—2020年長(zhǎng)江口徐六涇站年最高水位呈現(xiàn)出明顯下降趨勢(shì)，下降幅度為11.6mm/a，長(zhǎng)江口風(fēng)暴潮強(qiáng)度的減退是其主導(dǎo)因素。

（2）長(zhǎng)江口地區(qū)在臺(tái)風(fēng)登陸和發(fā)生極端風(fēng)暴高潮位的同時(shí)，遭遇極端降雨的概率較高，兩者同時(shí)超越各自數(shù)閾值（取95th百分位）的復(fù)合洪水事件重現(xiàn)期為89a。大通徑流與臺(tái)風(fēng)及其伴生事件在統(tǒng)計(jì)上近似獨(dú)立。

（3）忽略變量之間的相關(guān)性會(huì)造成復(fù)合洪水災(zāi)害風(fēng)險(xiǎn)的低估，且災(zāi)害因子重現(xiàn)期越長(zhǎng)，復(fù)合洪水災(zāi)害重現(xiàn)期低估程度越高。若不考慮相關(guān)性，風(fēng)暴潮“二碰頭”事件風(fēng)險(xiǎn)將被低估77.7%。