基于經驗模態分解重構方法的北江飛來峽洪水頻率分析

關鍵詞：變化環境；一致性；水文頻率；設計洪水；經驗模態分解；EMD；北江；飛來峽

中圖分類號：TV122 文獻標識碼：A 文章編號：1001-9235（2024）12-0001-10

水文頻率（概率）分布估計理論和方法是工程水文計算的核心科學問題［1］。傳統的水文頻率分析方法的一個基本前提是水文系列滿足一致性假設［2］，即水文極值的概率分布或統計規律在過去、現在和未來保持不變［3］。然而，美國水文學者Milly等［4］在Science雜志上提出了“一致性已經消亡（stationaryisdead）”的論斷，即在氣候變化和人類活動的影響下，水文序列已經不再滿足一致性［5］。IPCCAR6［6-7］指出，隨全球增暖幅度增加，強降水事件的頻次加速增長，并且越極端的強降水事件，其發生頻率的增長率就越大。比如，在全球2℃溫升下，當前20年一遇的強降水事件發生頻率將增加22%，100年一遇的強降水事件發生頻率將增加45%以上。采用傳統的工程水文方法制定的暴雨與洪水設計值，將面臨由變化環境帶來的風險［8］。

變化環境下水文頻率分析的核心內容包括水文序列非一致性檢驗方法與非一致性水文頻率分析方法2個方面［1-2］。在針對趨勢、跳躍等的非一致性檢驗方法中，Mann-Kendall、Hurst指數法和貝葉斯方法等非參數方法應用較為成功［2］。謝平等［9］提出了水文變異診斷系統，對各種方法的結果進行綜合，以期使變異檢驗結果更加可靠。非一致洪水頻率分析方法主要有：還原還現法、條件概率分布法、混合分布法和時變矩法等［10］。其中，還原/還現法應用最為廣泛，這種方法認為非一致性水文序列由確定性成分和隨機性成分構成，確定性成分（趨勢、跳躍等）通常被定義為非一致性成分；而隨機性成分定義為一致性成分，可通過傳統的頻率計算方法確定［11-13］。還原/還現法中確定性成分的分解和擬合是本方法的關鍵所在，通常使用最小二乘法等方法來擬合趨勢線。經驗模態分解方法被認為是200a來以傅立葉變換為基礎的線性和穩態頻譜分析的一個重大突破，在處理非線性非平穩信號方面具有非常明顯的優勢［12-13］。馮平等［14］、謝平等［15］、趙雪花等［16］引入經驗模態分解方法來對非一致性水文序列進行分解重構與頻率分析，表現出廣闊的應用前景。為此，本文采用了傳統的頻率分析方法和基于EMD（Empirical Mode Decomposition）分解重構方法對北江飛來峽站的洪水頻率進行了計算和比較分析，為確定北江洪水重現期與洪水設防標準提供參考。

1研究區及數據

北江是廣東省境內流域面積最大的江河，在思賢滘以上集水面積46710km2，其中廣東省境內42930km2，約占廣東省國土面積的1/4。北江干流至三水區思賢滘全長468m，其中，中游段為韶關市沙洲尾至清遠市飛來峽，河長173km，河道平均下游段坡降0.125‰，河谷多呈U字形；下游段為飛來峽至三水區思賢滘，河長83km，河道平均坡降0.0815‰，處于平原區，河面寬闊，兩岸多堤防［17］。北江流域的地勢北高南低，正處在西南季風的迎風坡，加之河流水系呈闊葉脈狀分布，洪水匯流集中迅猛，暴漲暴落［18-19］。北江的暴雨和洪水多發生在前汛期（4—6月），尤其以5、6月最為集中。北江流域平均降水量1807mm（1956—2016年），清遠、珠坑、橫石一帶是全省的降水高值區之一［21-20］。北江洪水主要來自橫石以上地區，下游防洪控制斷面石角站年最大洪水的15d洪量中，橫石站來量占84％。由于流域面積不大，一次較大的降雨過程幾乎可以籠罩整個流域，加之流域坡降較陡，橫石以上的干、支流洪水常常遭遇。橫石以下支流的發洪時間一般稍早于干流，較少與干流洪水遭遇［21］。最近短短31a（1994—2024年），北江相繼發生“94·6”（超50年一遇，按原洪水設計標準計，下同）、“97·7”（20年一遇）、“06·7”（50年一遇）、“13·8”（接近50年一遇）、“22·6”（超100年一遇）、“24·4”（接近100年一遇）等大洪水及特大洪水。這一方面說明北江流域的暴雨洪水特性可能發生了變異，極端洪水事件呈現出趨多趨頻趨廣趨強的態勢；另一方面也對原有的洪水頻率計算方法和設計洪水標準提出了質疑。

橫石水文站和飛來峽水文站都是國家基本水文站，分別位于飛來峽水利樞紐壩址的上游和下游，兩站相距7km，距壩址分別約5、2km，兩站集水面積相差約0.6%，區間無支流匯入，由于飛來峽水利樞紐的興建，橫石水文站于1999年1月遷移至飛來峽水文站，兩站的基本情況見表1，位置見圖1。

本文主要研究的是年最大洪峰流量，因橫石水文站和飛來峽水文站集水面積差異不大，視為一個站，統稱為飛來峽站，兩站的資料系列以1999年元旦為界前后相接。按照《珠江洪水調度方案》（2014年），飛來峽水庫當入庫流量大于5000m³/s但不大于15000m³/s（約20年一遇）時，水庫敞泄；當入庫流量大于15000m³，水庫控制出庫流量，起到了一定的滯洪調峰作用。可見，飛來峽水庫對于5000～15000m³/s的流量過程影響不大；而對于最高流量大于15000m³/s的特大洪水，1999年后共發生4次（“06·7”“13·8”“22·6”“24·4”），但由于飛來峽庫區英德防洪片為常住人口近30萬人的英德市城區，庫區內波羅坑、連江口、社崗防護片均存在淹沒損失大、啟用決策難的問題，制約著飛來峽水庫高水位調度運用［22］，因而飛來峽水庫發揮的滯洪調峰作用有限，飛來峽入庫流量與下游飛來峽站的實測流量相差不大。同時，北江流域除飛來峽水庫外，還存在著其他大量的水庫工程聯合調度，因此還原計算非常復雜。為此，對于“06·7”“13·8”“22·6”等3場特大洪水，本文直接采用了飛來峽站整編資料的實測洪峰流量；對于2024年剛發生的“24·4”特大洪水，由資料暫未整編，本文采用報汛資料中飛來峽水庫的入庫流量。所有的資料均來自于廣東省水文局。

飛來峽站（橫石站）的設計洪水值由水利部珠江水利委員會1983年核定并于2005年復核（表2），2次計算所采用的資料系列分別至1982、1997年，采用的P-Ⅲ型曲線的期望值、Cv和Cs分別為9450、0.34和1.02，2次的計算結果基本相同，本文稱之為“2005年頻率曲線”。

2研究方法

2.1水文變異診斷方法

使用曼肯德爾法（Mann-Kendall，M-K）來檢測水文序列的變化趨勢，并用重標極差分析法來檢測變化趨勢的持續性。曼肯德爾法是由曼（H.B.Mann）和肯德爾（M.G.Kendall）提出的一種用于檢測序列變化趨勢的非參數檢測方法［23］。這一方法的優點是不需要樣本遵從一定的分布，也不受少數異常值的干擾，適用于類型變量和順序變量，不僅可檢測序列的變化趨勢，還可以明確突變發生的時間及突變區域［24］。M-K檢驗的計算方法可參考相關文獻。

重標極差分析法（Rescaled Range Analysis，R/S）是英國水文學專家Hurst提出的一種非線性科學預測方法［25］，用于判別序列對時間的依賴性。根據統計量H∈［0，1］判斷序列的持續性：當H=0.5時，序列是標準的隨機游走，即現在對未來不會產生影響；當Hgt;0.5時，序列變化前后正相關，未來的趨勢與過去一致，H越接近1，持續性越強；Hlt;0.5時，序列變化前后負相關，未來的趨勢與過去相反，H越接近0，反持續性越強［26］。R/S檢驗的計算方法可參考相關文獻。

2.2水文頻率計算方法

洪水資料的選樣采用年最大值法，即每年選取一個最大洪峰流量。洪水頻率曲線線型選用了皮爾遜Ⅲ型（P-Ⅲ型）曲線。P-Ⅲ型曲線有期望值Ex、變差系數Cv和偏態系數Cs等3個參數，本文中的Cs/Cv沿用了以往本站頻率計算中采用的經驗值3.0。使用矩法初估參數，通過計算機程序按離差平方和最小二乘準則（Ordinary Least Squares，OLS）優化適線，得到最優參數估計值，在此基礎上計算得到各頻率的洪水設計值［27-28］。經驗頻率的計算采用統一處理法［29］，將實測洪水和歷史洪水共同組成一個不連序的系列，認為它們共同組成一個歷史調查期為N年的樣本，在N年中統一排序。

2.3經驗模態分解方法

經驗模態分解方法（EMD）是美國NASA的Huang等［12］于1998年創造性地提出的一種新型自適應信號時頻處理方法，特別適用于非線性非平穩信號的分析處理［30］。EMD的本質是將原始時間序列分解為若干IMF（Instinct Mode Function，本征模態函數）分量以及殘差項（Residual Error，RES），即由式（1）組成：

式中：x（n）為原始時間序列；n為時間序列的長度；IMFi（n）為各個本征模態函數分量；Res（n）為殘差。IMF分量是具有時變頻率的震蕩函數，IMF分量的震蕩頻率總體上隨著階數i的增大由高頻變為低頻，分別代表一組特征尺度（頻率）的數據序列，它反映非平穩信號的局部特征。RES又稱為趨勢項，反映非平穩信號的整體變化趨勢。任何一個時間序列或信號都由長期趨勢、各種周期性成分以及隨機成分組成，水文氣象要素時間序列也不例外。

通過EMD分解，周期性成分與隨機成分被分解到各階IMF分量中，可用傳統的頻率分析方法計算重現期；而長期趨勢成分則被分離出來成為RES項，可通過擬合方程來分析其發展變化。EMD分解的具體計算步驟方法可參照文獻［31］，本文借助Matlab軟件的EMD函數進行了分解。

3結果與討論

3.1實測歷年最大洪峰流量變化過程及趨勢性檢驗

飛來峽站1953—2024年（其中，1953—1998年為橫石站，下同）72a間歷年實測最大洪峰流量過程線見圖2。根據水利部《全國主要江河洪水編號規定》（2019年），當北江石角水文站流量達到12000m3/s時，進行北江洪水編號，本文參照這一流量值對飛來峽站的洪峰流量進行評價。1953—2024年72a間飛來峽站發生編號洪水的年份有16a。其中，1953—1993年，41a的時間里飛來峽站發生編號洪水的年份共有5a，分別是1955、1964、1968、1972和1983年，平均每8a才有1a會發生編號洪水；而1994—2024年，31a的時間里編號洪水的年份達到了11個，平均不到3a就有1a會發生編號洪水。20a滑動平均線的數值在1990年代及之前不足10000m³/s，在2000年前后突破10000m³/s，最近10a更是突破了12000m³/s，前后上升幅度達到了3000m³/s。由此可見，飛來峽站極端洪峰流量呈現出趨頻趨大趨多的特征。

利用曼肯德爾法（M-K）對歷年最大洪峰流量進行趨勢性檢驗和突變性檢驗。M-K趨勢性檢驗的統計量Zc為2.35，大于正態分布顯著性水平P=0.05對應的臨界值（1.96），表明72a間年最大洪峰流量總體呈顯著上升趨勢。

M-K突變性結果見圖3。UFk線和UBk線有2次交叉，第一次交叉于1984年，UFk線向下穿過UBk線后繼續向下發展，直到1991年才重拾升勢，表明這一時期的年最大洪峰流量呈現下降趨勢，這與圖2中1984—1991年8a間所有的年最大洪峰流量不到10000m³/s的現象是一致的；第二次交叉于1994年，UFk向上穿過UBk線，并從2004年開始UFk呈快速上升趨勢，在2013年后向上突破1.0后繼續波動上升，在2022年向上突破了顯著性水平P=0.05的臨界值（1.96），表明飛來峽的年最大洪峰流量在近30a間發生了突變，突變的起點為1994年（“94·6”洪水）。

利用R/S方法計算飛來峽站的歷年最大洪峰流量系列的Hurst指數為0.7141，遠大于0.5，表示時間序列具有長期記憶性，即未來的變化趨勢將繼續保持向上。

3.2飛來峽站歷史洪水的量級與重現期考證

水利部珠江水利委員會于1983年協同廣東省三防指揮部辦公室、廣東省水文總站等單位對北江橫石站的歷史大洪水進行了考證和定級［32-33］。從1983—2024年已有41a，飛來峽站的實測資料系列（1953—2024）比1982年（1953—1982）時延長了1倍有余，期間又發生了多次與考證洪水量級相當的特大洪水和大洪水。本文根據前述的考證資料和自1953年設站以來的實測資料，確定飛來峽站歷史洪水排位和重現期見圖4，主要分析結果如下：①1915年考證洪峰流量為21000m³/s，為1764年以來首位，重現期為261a（1764—2024年）；②2022年實測洪峰流量為20600m³/s，為1764年以來第2位，重現期為130a；③1931年考證洪峰流量19600m³/s和2024年入庫流量18900m³/s，分列第3、4位，重現期分別為87、65a；④1982、1764、1877年洪峰流量同等量級，約為18000m³/s，排在1764年以來第5—7位，重現期約為50a；⑤1878、1914年考證洪峰流量低于1764年，但未有具體數值，可能僅與1994年相當，重現期低于40a；⑥1994年洪峰流量約為17500m³/s，排在歷史洪水第8位或第8位之后，排在實測洪水第4位，重現期約為35a，按常遇洪水處理；⑦2006、1997、2013、1968年的洪峰流量分別為16300、15900、15800、15000m³/s，排在實測洪水第5—8位，按常遇洪水處理，重現期約為20～25a。

與前述1983年的考證排序相比，本次排序主要變化是增加了排在第2位的“22·6”洪水和排在第4位的“24·4”洪水，并將序列長度從219a延長到261a，使得原有洪水的名次和重現期發生了變化。其中：1915年大洪水的重現期因此由219a增加至261a；1931年洪水的重現期由第2位下移至第3位，重現期由110a降低為87a；“82·5”洪水，1877、1764、1878、1914年等洪水排名退后2位，重現期相應變小。

3.3考慮歷史洪水的頻率計算

使用統一處理法進行經驗頻率計算。歷史調查期為1764—2024年共261a；發生歷史特大洪水的年數共7a，按洪峰流量從大到小排列為1915、2022、1931、2024、1982、1877、1764年；實測期年數為1953—2024年共72a；其中發生在實測期的特大洪水的年數有3a，分別為2022、2024和1982年。前7位洪水按特大洪水計算經驗頻率；其余68a的洪水按一般洪水計算經驗頻率。計算結果見表3。

經計算得到，飛來峽站歷年洪峰流量的均值Ex的初步估計值為10029m³/s，Cv的初步估計值為0.33。第3個參數偏態系數Cs采用了北江流域的經驗取值3.0Cv即0.99。經適線并按OLS原則優選參數，得到Ex、Cv和Cs的最優估計值分別是10083m³/s、0.32和0.96，各頻率的洪水設計值見表4，擬合曲線見圖5。

與2005年頻率曲線相比，均值由9450變為10083m³/s，增大6.7%；Cv值由0.34變為0.32，曲線整體向上移動，各頻率設計值增大500～600m³/s。

3.4 基于EMD分解與合成的頻率計算方法

采用Matlab軟件EMD函數對1953—2024年最大洪峰系列進行分解，得到4個IMF分量與1個趨勢項RES，見圖6。由于歷史洪水不連序，在本方法中不予考慮。從圖6中可以看出：各IMF分量圍繞0軸上下波動，頻率由高到低依次遞減，幅值變化也是依次遞減；而趨勢項RES呈單調上升趨勢，并在1990年代開始上升速率加快，變化區間在9000～15000m³/s。

將各IMF分量進行疊加，記作IMFz，IMFz圍繞橫向0軸上下波動，無明顯的趨勢性，即具有較好的一致性，見圖7a。由于P-Ⅲ型頻率分布曲線須保證數據序列為非負數，因而不能直接使用IMFz序列進行分析。為此，將IMF_z序列與原序列的均值相加得到一新序列，記作（fx），用以表征原序列的一致性成分，見圖7c。將RES項進行曲線擬合，選用四次多項式取得了非常好的擬合效果，效率系數達到0.999，見圖7b，并且曲線表示的是年最大洪峰流量整體變化特征，取值全部為正值且高懸于橫向0軸之上。為此，將RES擬合曲線方程減去原序列的均值，得到一個新的序列g（x），用以表征原序列的趨勢性成分，見圖7d。

這樣，原始序列分解成隨機性成分f（x）和趨勢性成分g（x）2個子序列。將f（x）子序列擬合P-Ⅲ型頻率分布曲線，估計得到的擬合參數：E_x為9485，C_v為0.35，C_s為C_v的3倍即1.05。此頻率分布曲線各頻率對應的年最大洪峰流量計算結果見表5中“過去平均”一行。根據趨勢性方程g（x）計算各代表性水平年的趨勢性成分，1960、1980、2000、2020、2023年分別為-700、-400、1300、3500、3800m3/s

3.5 分析討論

北江飛來峽站年最大洪峰流量2005年頻率曲線、2024年頻率曲線和采用基于EMD分解重構方法的頻率曲線計算參數見表6，頻率曲線見圖8，其中EMD方法只列了1980、2000、2020年等3個水平年。表和圖中可以看出，五者之間Cv的變化范圍為0.32～0.35，變化幅度較小；主要區別在于均值Ex及趨勢性成分的變化。2024年頻率曲線（即設計值）比2005年頻率曲線整體上移約600m³/s；EMD1980與2005年頻率曲線在小頻率時（Plt;5%）時基本重合，在大頻率時（Pgt;5%）EMD1980年頻率曲線略低；而EMD2000、2020年曲線分別比2005年頻率曲線整體上移約1300、3500m³/s。當P=0.5%（200年一遇）時，5條曲線對應的設計值分別為21000、21200、20600、22300、24500m³/s，后兩者要遠大于前三者。而對于同一量級的洪水，如21000m³/s（1915年量級），對應5條頻率曲線的重現期分別為200a、約200a、約200a、略超100a和接近50a。也就是說，采用傳統方法的2024年頻率曲線與2005年相比變化不大，而采用EMD分解重構方法的頻率曲線在1980年水平年時的變化尚不顯著，但在2000、2020年水平年卻發生了顯著變化，直接影響到了洪水的定級。再對照圖2中飛來峽站年最大洪峰流量的20a滑動平均線可知，2000年前后的滑動平均平均值比1990年前上升1000m³/s以上，而2020年前后的滑動平均值比1990年前上升3000m³/s以上，這與EMD方法計算得到的趨勢性成分基本一致。由此可見，EMD分解重構頻率計算方法能夠及時地反映出變化環境下的趨勢性成分變化特征，同時能夠較好地解釋關于近年來頻繁出現的多次極端洪水的疑問，具有一定的合理性和推廣應用價值。

傳統的洪水頻率計算方法（如2005、2024年頻率曲線）的一個重要前提是水文事件的一致性假設，即假設過去、現在及未來的水文樣本來自同一總體。但是，在氣候變化和人類活動的雙重影響下，部分地區的水文規律已發生顯著變化，并且年代相隔越久遠，變化幅度可能越大。因此，部分地區的水文一致性假設已不再成立，迫切需要探索使用變化環境下水文頻率計算方法，其中，本文使用的基于EMD重構的水文頻率計算方法是一個值得深入研究的方向。同時，無論是使用傳統的頻率分析方法，還是使用基于EMD的水文頻率分析方法，北江飛來峽站的水文頻率計算結果都發生了較大的改變，因而需要盡快重新核定北江的設計洪水數值，確保防洪工程體系和調度決策體系能夠具備足夠的安全保障能力。

4 結論

a）北江極端洪水自1990年代中期表現出趨多趨頻趨強的特點，經檢驗發生了顯著變異，并且這種變異具有一定的持續性。

b）洪水資料系列的延長與近年來發生的多場大洪水，改變了北江洪水的排序與重現期。按傳統的水文頻率計算方法重新適線計算后，頻率曲線發生顯著上移，各頻率洪水設計值需比原設計值增大500～600m³/s。

c）使用基于EMD分解與合成的頻率計算方法，在考慮趨勢成分后，2020水平年的各頻率設計值比原設計值增大3500～3800m³/s，2030水平年的各頻率設計值比原設計值增大3800m³/s以上。北江“22·6”和“24·4”等原超100年甚至超200年一遇的洪水變成不到50年一遇。

d）基于EMD分解合成的頻率計算方法更能反映出變化環境下洪水頻率的變化特征，具有一定的合理性和推廣應用價值。

e）北江飛來峽站的水文頻率計算結果都發生了較大的改變，需盡快重新核定。

f）EMD方法的局限性主要包括模態混疊和末端效應，其中末端效應可能會使本文的EMD分解重構方法會產生不合理結果，即在數據序列的開始和結束部分，由于沒有足夠的數據點進行考慮，極值包絡線在端點處會發生發散，產生較大的誤差，為此，可通過EEMD（Ensemble Empirical Mode Decomposition，集合經驗模態分解）分解、CEEMD（Complementary Ensemble Empirical Mode Decomposition，互補集合經驗模態分解）等方法加以改進。同時，在趨勢分析中，現有的趨勢在未來時段能否繼續保持，也是一個不確定的問題，因此，本文所介紹的方法趨勢外延的時段不宜過長。