珠江河口區復合洪水事件多要素非一致性的時變遭遇分析

張 卡，劉丙軍,2，胡仕焜，曾 慧，張明珠，李 丹

（1.中山大學 土木工程學院，廣東 珠海 519085；2.中山大學 水資源與環境研究中心，廣州 510275；3.廣州市水務科學研究院有限公司，廣州 510220）

近年來，中國南海沿岸極端臺風發生的頻率與強度逐年增大，上游洪水、臺風和天文潮疊加引起的復合洪水事件問題日趨嚴重，嚴重威脅三角洲地區人民的生命財產安全（羅志發 等，2022）。上游洪水、臺風和天文潮的時空分布不均，引發的潮位增水，洪水、狂風等事件組合存在高度不確定性。針對自然事件的風險特征，開展上游洪水、臺風和天文潮復合事件的主控因子診斷及其和多要素遭遇分析，對河口地區防范復合洪水事件有重大的現實意義。

復合洪水事件多要素遭遇分析是研究復合洪水事件規律的重要手段之一，利用Copula函數連接多個變量邊緣分布構建聯合分布函數的水文統計方法，既能通過篩選不同的連接函數體現隨機變量之間關聯性，又能通過構造多變量不同的邊緣分布（關帥 等，2015），使得每種邊緣分布函數保留自己分布的特點，反映單個隨機變量分布特征。Copula函數可以反映水文事件隨機性、相依性的特點，近年來被廣泛應用到水文要素的頻率分析（趙鐵松等，2020）。目前，國內基于上海地區降水、風速、潮位、雨量的年極值時間序列，探討了“雨洪風”“雨洪潮”和“雨風潮”3種致災因子碰頭情況下的聯合分布函數，并進行了組合風險分析（賀芳芳等，2021）。高月嬌等（2023）構建降水廣義可加模型的二維Copula模型分析得出，1982—2015年黃土高原地區旱澇復合事件極端情景顯著上升，且復合事件動態變化的主導因子為北極濤動指數和太陽黑子指數；侯靜惟等（2019）構建了熱帶氣旋災害的主要致災因子大風和降水之間的二元Copula 函數，研究了不同Copula的聯合重現期的相關性。許瀚卿等（2022）通過構建上海極端降水和極值水位兩變量的聯合分布函數，評估不同聯合重現期下降水和潮汐增水組合的復合洪澇風險。Bevacqua 等（2019）利用Copula 結構研究發現，意大利拉文納2015年2月的復合洪水是通過一個低壓系統在多個河流集水區產生臺風、天文潮和強降水形成的。Ohlwein 和Friederichs（2008）基于Copula 的方法運用多元關聯分析工具箱（MvCAT）、Vine Copula中C藤和D藤、貝葉斯概率等方法，擴展了多元模型的數量。Olbert 等（2023）以愛爾蘭南海岸的科克市為例，將統計和流體動力學模型聯系起來，利用Copula函數確定浪涌、沿海水位、河流流量關系的聯合概率，提高了MSN_Flood模型模擬不同洪水情景下洪水傳播情況的精度。Yang 和Qian 等（2019）提出了一種有效的粒子群優化（PSO）估計風速、極端水位、降水邊際累積分布和Copula函數的參數，該方法計算Copula 條件聯合分布誤差較小。目前利用二元Copula函數對復合事件進行遭遇分析的方法成熟且研究廣泛，這些研究表明極端復合事件發生概率正顯著上升，復合事件中主要致災因子的影響程度不同。受氣候變暖的影響，復合事件中各個致災因子并非呈現出一致性，且復合事件的發生受多種致災因子影響，二元Copula函數方法應用受限，未來仍需深入研究多元Copula函數技術以及復合事件中非一致性水文資料的遭遇分析。

當前對沿海地區復合洪水事件遭遇分析的研究大多以水文因子一致性假設為基礎，要求水文序列的統計特征不隨時間改變。珠江河口河網縱橫交錯，天文大潮、臺風等極端氣候背景和人類活動雙驅動作用下，會導致對洪水概率分布描述的偏差（顧西輝 等，2014）。傳統以固定Copula 參數分布的研究方法，忽略了各個水文序列的非一致性，會導致復合洪水事件遭遇等級有所偏差。因此，本研究以珠江河口為研究對象，在傳統復合洪水事件聯合分布風險評估框架內，進一步考慮潮位、風速、流量3 種因子潛在的非一致性，并利用GAMLSS模型（Durocher et al., 2015） 和Copula 似然比方法（Huang et al.,2017）構建三變量聯合時變分布函數，分析不同等級復合洪水事件的時變特征和遭遇規律。以期為珠江河口復合洪水事件防御提供理論與技術支撐。

1 研究區概況與數據

珠江河口位于廣東省中南部、珠江下游，該地區總面積為5.54 萬km2，珠江干流全長2 214 km，主要由西江、北江和東江匯合而成。水道縱橫交錯形成密集的河網流經虎門、蕉門、洪奇門、橫門、磨刀門、雞鳴門、虎跳門和崖口八大口門匯入南海。這些密集的河網攜帶的泥沙在珠江口河口灣內堆積形成了包括惠州市、東莞市、廣州市、佛山市、江門市、中山市、深圳市、珠海市、肇慶市九大城市的珠江三角洲（圖1）。珠江入海的水量約為3 742 億m3，其中，西江是珠江的主干流，西江、北江的中下游干流水道主要從珠江的磨刀門入海，三水站和馬口站分別位于北江和西江中下游干流處，可代表珠江上游總徑流量，而澳門站和三灶站位于珠江河口，能較好的代表珠江河口地區受臺風影響時的潮位和風速。每年汛期西北太平洋海面發生的臺風大多在珠江河口地區登陸，若同時和天文潮、上游洪水疊加往往會造成嚴重的水位抬升，導致漫灘漫堤，造成重大的生命財產損失（李闊 等，2010）。2018-09-07“山竹”臺風在西太平洋海面生成，9月16日在中國廣東省臺山市海晏鎮登陸，引發珠江河口測站最大增水達3.37 m，造成廣東省5人死亡，約252 萬人緊急避險轉移和安置，600 余間房屋損毀，103.5千hm2農作物受災，直接經濟損失32.1 億元（劉士誠 等，2021）。粵港澳大灣區位于珠江河口地區，地理位置特殊，經濟發達，受上游洪水、臺風和天文潮復合洪水事件影響巨大。

圖1 研究區域示意Fig.1 Schematic diagram of study area

基于珠江河口地帶1988—2017年共30 a的三灶站每日最大潮位、澳門氣象站每日平均風速，采用年最大值法篩選三灶站日極值潮位，并找到該天文潮事件對應的日最大風速、日最大流量。其中，極值潮位表征珠江河口地區復合洪水事件的天文潮量級，三水和馬口兩站每日總流量代表復合洪水事件的洪水洪峰量，澳門日最大風速表征復合洪水事件的臺風量級構建三要素時變聯合分布函數，分析珠江河口3種因子時變特征及其遭遇規律。

2 研究方法

2.1 GAMLSS函數

以時變矩法的內容作為基本思路，Rigby 等（2005）提出水文序列中位置、尺度和形狀的廣義可加模型（Generalized Additive Models for Location, Scale and Shape, GAMLSS），用于模擬不同水文序列的時變參數邊緣分布。模型定義研究變量Y的n個獨立觀測值yi服從概率分布f(yi|Θi)，其中Θi=(θi1，θi2，θi3，θi4)=(μi，σi，νi，τi)為概率分布參數組成的向量，μi和σi對應為變量Y的均值向量以及變差系數向量，νi和τi定義為偏度和峰度參數。利用連接函數gk(θk)建立概率分布參數和解釋變量之間的聯系，連接函數的一般表達式為（Rigby et al., 2005）：

式中：gk(θk)代表n個時段的解釋變量Xk對已知設計矩陣Zjkτμ的連接函數，k= 1，2，3，4；βk為待求解的參數向量；Zjk為已知的設計矩陣；τμ是由隨機變量組成的列向量，選取不同的概率分布曲線和概率分布參數可以構建多個GAMLSS 模型，本文GAMLSS 模型主要基于R 語言中的GAMLSS 包，通過其中可選擇的伽馬分布（Gamma, GA）、耿貝爾分布（Gumbel, GU）、邏輯斯諦分布（Logistic,LO）、 正 態 分 布（Normal, NO）、 韋 伯 分 布（Weibull, WEI）、廣義帕累托分布（Generalized pareto, GP）等選取最優時變邊緣分布曲線。

2.2 時變Copula函數

不同類型Copula函數與參數分別代表多個變量耦合的形狀以及耦合強度?；贑opula函數的似然比方法（Copula-based likelihood-rat io test, CLR test）可以檢驗不同變量耦合強度（Huang et al.,2017）。假設Copula 函數的類型不變，檢測Copula函數的參數是否發生變化，其原理為：已知yi是序列觀測值，且服從概率分布函數F(y1|θc，n)，其中θc,n為不同時段的Copula函數的參數，H0假設Copula函數不同時段的參數不發生變化。

構造與Copula函數參數有關的似然比Xk，當Xk較小時，拒絕原假設H0，兩者擬合而成的Copula函數存在變點K， 似然比Xk與Copula 函數參數θc，0，θc，1，θc，2極大似然值的解析式為（方偉，2020）：

式中：Ln()為數據數量為n的水文序列似然函數；Lk()和Ln-k()為Copula 突變點前后的似然函數；c()為Copula 函數的概率密度函數；η0，η1，η2為Copula函數參數θc,0，θc,1，θc,2的極大似然估計值。

若Xk的對數統計檢驗值Zn概率值＜5%，表明在5%顯著性水平下拒絕原假設，即該變量的Copula關系存在突變點，其突變點位于（方偉，2020）：

若Copula 函數存在突變則其參數θk為時變參數。當Copula 參數隨時間發生變化時，引入年份t作為解釋變量，其解析式為（方偉，2020）：

采用邊際函數推斷法（Inference function for Marginal method, IFM），通過最大化對數似然公式求Copula參數θ的表達式（方偉，2020）。

式中：c為Copula的概率密度函數；f1、f2分別代表潮位、風速的已知邊緣分布函數，將等號右側第一項最大化即可得到Copula時變參數的最優β0和β1。

Joe（1993）引入嵌套阿基米德Copula 即分層Copula，來改善不同變量的相關性。本文應用完全嵌套Copula 函數（陳沖，2020），將潮位、風速二變量時變Copula函數與和時變流量邊緣分布函數進行擬合，構建時變三變量聯合分布函數，其表達式為（張野 等，2017）：

式中：un指不同變量的邊緣分布函數；C表示不同的Copula連接函數。

3 結果分析

3.1 時變聯合分布函數優選

根據單位根檢驗（Augmented Dickey-Fuller test）得到潮位、風速、流量時間序列均為非一致性序列，選用GAMLSS模型估計3組水文時間序列的時變統計參數，利用貝葉斯信息度量（Bayesian information criterion）BIC值，判斷不同邊緣分布曲線的擬合程度，如表1 所示。BIC 值最小的Normal分布、Gamma分布和Weibull分布的擬合程度越高，它們分別為潮位、風速、流量擬合的最優邊緣分布函數。圖2顯示：擬合后的邊緣分布函數偏差點均位于95%置信區間（虛線）內，且均靠近偏值為0的黑色虛線周圍，即3種變量非一致性邊緣分布函數擬合效果較好，擬合后的邊緣分布函數各年份偏差值較小，擬合誤差小，精度高，能較好地模擬不同水文序列概率密度的時變特征。

表1 GAMLSS模型擬合的潮位、風速、流量BIC值Table 1 Tide level, wind speed and flow BIC values fitted by GAMLSS model

圖2 基于Worm plot的潮位（a）、風速（b）、流量（c）擬合優度檢驗Fig.2 Tide(a), wind speed(b) and flow(c) fitting effect test based on Worm Plot

潮位和風速的相依關系呈現非一致性特征，利用Copula 似然比方法進行擬合。選取BIC 最小的Frank Copula 函數作為最優的Copula 函數。優選的Frank Copula函數參數θ和時間的函數關系為ln (θ-1)= 0.84 -0.54t（其中t為時間年份與初始時間1988年之間的差值），構建時變參數的潮位-風速時變Copula參數（圖3）。在潮位-風速時變Copula參數圖中，當潮位、風速耦合結構呈現非一致性關系時，選取的最優Copula的參數隨時間不斷減小，由于Copula 函數類型均為Frank Copula，潮位和風速相依結構的形狀不變，兩者尾部相關性不強。選取不同年份的潮位-風速聯合分布函數，與時變最優流量邊緣分布函數構建第三層Copula函數，通過嵌套結構獲得1990、1999、2017年的三變量聯合分布函數（圖4），可知，隨著時間變化3種變量耦合后的聯合分布的概率有小幅度增大。

圖3 潮位-風速時變Copula參數Fig.3 Time-varying Copula parameter diagram of combined tidal level and wind speed

圖4 各年三變量聯合分布函數（a.1990年，θ=3.17；b.1999年，θ=2.94；c.2017年，θ=2.54）Fig.4 Three-variable joint distribution function in 1990(a.θ=3.17), 1999(b.θ=2.94), 2017(c.θ=2.54)

3.2 不同等級臺風重現期時變分析

利用非線性Copula函數，分析2008-09-19—26“黑格比”、2012-07-21—24“韋森特”與2013-09-18—222“天兔”3場不同等級復合洪水事件在1988—2017年不同年份時變聯合重現期，利用三場等級復合洪水事件的潮位、流量、風速極值計算時變聯合重現期，結果如圖5所示，通過分析得出如下幾點認識：

圖5 不同臺風事件時變聯合重現期Fig.5 Combined return periods of different storm surge events

1）時變聯合分布函數能更好地體現水文序列的非穩定性和時變特征，能較好地反映隨時間變化下各等級臺風危險性的變化。假設水文序列為平穩序列，利用傳統的邊緣分布函數和傳統恒定Copula參數可以計算得到“黑格比”“韋森特”“天兔”復合洪水事件的恒定參數聯合重現期分別為11.9、4.37、6.51 a。但利用Copula 似然比方法計算的時變情況下1988—2017年“黑格比”復合洪水事件平均約為12.28 a一遇，2008年“黑格比”復合洪水事件的重現期為11.74 a，比傳統恒定Copula參數計算的聯合重現期小；利用恒定參數的聯合分布函數求得的“韋森特”事件“天兔”事件的重現期分別為4.37 和6.51 a，比復合洪水事件發生年份的時變聯合分布函數的重現期4.19 和6.53 a 長，且3 場復合洪水事件重現期隨時間變化不斷縮小，利用時變聯合分布函數比傳統恒定Copula參數方法計算的3場典型復合事件的重現期更長，即復合洪水事件等級較強，利用時變聯合分布函數比傳統恒定Copula參數計算的重現期更符合實際。

2）在變化環境下，隨著時間推移，不同等級復合洪水事件聯合重現期不斷縮短，發生的頻率不斷增大，且等級較大的復合洪水事件重現期減小的趨勢更快，說明未來相同復合事件發生頻率將會不斷增加，同一場復合洪水事件在未來發生的危險性增強。1999年發生“黑格比”等級復合洪水事件的重現期約為13.06 a，2017年該等級復合洪水事件的重現期為10.80 a，且在1999—2004 年其復合洪水事件重現期縮短較快，表明該等級復合洪水事件未來發生的頻率在增大，未來大約每12 a將可能出現一場該等級的復合洪水事件。“黑格比”“韋森特”“天兔”3 場復合洪水事件重現期在1988—2017 年的變化特征一致，在1988—2017 年“黑格比”“韋森特”天兔“復合洪水事件的重現期均縮短了2.26、1.06、0.97 a，整體上，等級越高的復合洪水事件縮短的趨勢越快，則未來相同等級的復合洪水事件重現期將會縮短，發生的概率將會增大。

3.3 復合洪水事件主控因子分析

基于珠江河口潮位、風速、流量的非平穩性變化特點與復合洪水事件遭遇組合，利用敏感性分析方法，探討該河口區不同重現期基準下上游洪水、臺風和天文潮復合洪水事件的主控因子，若潮位、風速、流量單變量重現期均為5、10、20、50 a 一遇時，控制潮位、風速、流量中2 種水文變量重現期不變，將某一因子的重現期提高到10、20、50、100、200、500、1 000 a 一遇，計算該致災因子重現期提高后的2013 年復合洪水事件的同現重現期（圖6），通過分析4幅熱力圖得出如下幾點認識：

圖6 以5、10、20、50 a為基準單變量不同頻率下三變量同現重現期Fig.6 The recurrence period of three variables with one variable and different frequencies by 5, 10, 20, 50 a

1）當以5 a為基準重現期時，復合洪水事件主要受臺風影響較大；以10、20、50 a為基準重現期時，重現期逐漸增大，復合洪水事件主要受天文潮影響較大。整體上天文潮對復合洪水事件的影響程度最大。潮位和流量2種因子的重現期不變均為5 a一遇時，降低風速因子發生頻率，使其風速單變量重現期增長，則復合洪水事件危險性變化最強。風速重現期從5 a一遇增大到100 a一遇時，復合洪水事件的同現重現期從11.91 a 一遇增大到42.94 a 一遇，其增大幅度為260.54%，表明復合洪水事件對天文潮的敏感程度更高。當潮位的重現期從5 a 一遇提高到100 a 一遇時，復合洪水事件同現重現期從11.91 a 一遇變為36.47 a 一遇，整體變化幅度為206.21%；但隨著上游流量的發生頻率減小，流量重現期增大，復合洪水事件同現重現期的變化程度較小，當潮位和風速為5 a一遇，流量為100 a一遇時，復合洪水事件的同現重現期僅為18.59 a一遇，整體變化幅度為56.09%。當以10、20、50 a為基準重現期時，潮位重現期改變為1 000 a一遇時，三變量同現重現期為323.72、357.54、494.38 a，比流量和風速改變導致的同現重現期更長，表明復合洪水事件對天文潮更為敏感，更容易受其影響。

2）重現期基準改變時其同現重現期改變幅度也有所不同，當重現期基準為5 a時，風速從5 a一遇改變為1 000 a 一遇時，同現重現期從11.91 a 增長到81.28 a，重現期增長了69.37 a；重現期基準為50 a時，隨著潮位這一單變量重現期增長，同現重現期增長到420.17 a，增長幅度最大。因此當3種因子的基準較大時，某一變量重現期改變，三變量同現重現期變化幅度最大。

4 結論

基于珠江河口1988—2017年三灶站潮位、澳門氣象站風速、三水和馬口流量的逐日極值數據，利用GAMLSS 模型和時變Copula 函數構建了復合洪水事件發生時3種因子的時變聯合分布函數，分析了變化環境下上游洪水、臺風和天文潮3種事件的遭遇分析以及該復合事件的主控因子，得到以下主要結論：

1）復合洪水事件中潮位、風速、流量水文序列呈現非平穩特征，且互相之間有一定的相依關系，因此，在進行復合洪水事件遭遇評估時，需考慮3種因子的非平穩性特征。本文所構建的時變聯合分布函數能較好地反映3種水文序列的非平穩特征，對研究復合洪水事件多要素的遭遇特征適用程度更高，更符合實際情況。

2）與傳統恒定Copula參數的聯合重現期相比，非線性時變聯合分布函數計算的3次典型復合洪水事件的重現期隨時間推移不斷減小，表明該等級復合洪水事件發生概率不斷增大，未來該等級復合洪水事件的周期性會不斷縮小，重現期將縮短，即未來復合洪水事件將增多。且等級越高的復合洪水事件重現期縮短的趨勢越快。如“黑格比”臺風引起的復合洪水事件在1988—2017年的聯合重現期隨著時間演變呈現不斷減小的趨勢，從13.06 a一遇，減小至10.80 a一遇，未來大約每12 a將可能出現一場該等級的復合洪水事件。

3）當3種因子重現期較小時，復合洪水事件主要受臺風影響較大；重現期逐漸增大，復合洪水事件主要受天文潮影響較大。當其余2種致災因子均為5 a 一遇時，隨著風速和潮位的發生等級增長到100 a一遇，其同現重現期增長到42.94和36.47 a一遇，增長了31.03 和21.56a，流量所引起的聯合重現期的增幅僅為56.09%，即天文潮和臺風對上游洪水、臺風和天文潮的復合洪水事件規模的影響程度較大，上游洪水的影響作用相對偏小。當以10、20、50 a為基準重現期時，隨著潮位重現期改變為1 000 a一遇，三變量同現重現期增長幅度最大，復合洪水事件更容易受天文潮影響。