基于完備集合經驗模態分解的贛江徑流多尺度變化特征

萬智巍，賈玉連，洪祎君，蔣梅鑫

(1.江西師范大學地理與環境學院鄱陽湖濕地與流域研究教育部重點實驗室，江西南昌330022；2.中國科學院地理科學與資源研究所陸地表層格局與模擬院重點實驗室，北京100101)

0 引 言

近年來隨著全球極端氣候水文事件的頻繁出現，深入探討氣候變化背景下的水文氣象響應過程成為相關國際研究計劃的重要方向之一[1]。這其中有關河川徑流長時間、多尺度變化周期和演變規律的研究成為人們關注的重點[2]。河流作為所在區域的集水、匯水的樞紐和中心，其徑流變化反映的是整個區域的水文變化與響應過程，因此成為全球變化及其影響的指示器[3]。

贛江[4-5]是鄱陽湖的主要供水來源，流域面積占到整個鄱陽湖流域的50%以上[6]。近年來隨著三峽水利工程的完工和蓄水，鄱陽湖補給長江明顯，鄱陽湖水位持續偏低[7]。相關研究表明，三峽水利工程運行后，在河道沖刷、來水減少和蓄水等因素的共同作用下，鄱陽湖的枯水季節相當于提前了1個月左右[8]。因此，研究贛江徑流變化的基本特征及其長時間變化趨勢和規律，不僅具有重要的科學意義，也是制定區域可持續發展戰略迫切需要解決的現實問題。

有關研究還表明，水文系統同時具有多層次性和多時頻、多域頻變化結構[9]。因此，利用傳統的統計分析方法難以揭示水文氣象序列中所隱含的規律性變化過程[10]。為此，除了利用更長時間的水文觀測記錄進行統計分析外，經驗模態分解(Empirical Mode Decomposition，EMD)方法被認為是一種可以分解非線性復雜信號序列的新技術。EMD方法基于Hilbert-Huang變換，利用各信號周期不同的振幅頻率特征，將原始序列分解為一系列具有不同時間尺度的信號，即本征模函數(Intrinsic Mode Function，IMF)。由于各IMF信號具有很高的信噪比，因此這一方法被廣泛應用于水文、氣象、地球物理等領域的多尺度時頻分析[11-13]。本文利用最新提出的EMD改進方法——完備集合經驗模態分解[14-16](Complete Ensemble Empirical Mode Decomposition，CEEMD)，對贛江有器測數據以來的徑流序列進行多尺度變化特征和規律研究。

1 數據與方法

1.1 徑流數據

贛江的徑流數據來自1950年～2015年南昌外洲水文站的實測數據。該水文站位于江西省南昌市境內的贛江下游右岸，集水面積為80 948 km2，占整個贛江流域的99.6%[17]。外洲站設立于1949年10月，自1950年開始有連續的水文觀測記錄，是目前贛江流域具有最長連續水文觀測資料的站點。

1.2 研究方法

利用EMD的最新改進的CEEMD方法提取徑流序列中的周期和趨勢項，確定變化隱含的周期和最終變化趨勢，具體流程為[14-16]：首先，原始信號中加入2k次均值為0，符號相反的常數白噪聲Ni(i=1，2，…，2k)，獲得新的時間序列fi(t)；其次，對fi(t)進行EMD分解，得到各IMF分量；最后，得到原始序列

(1)

式中，cij(t)為第i次加入白噪聲后得到的第j個分量，也稱之為IMF(本征模函數)；rij(t)為第j次分解信號后得到的剩余趨勢項。

利用線性傾向和MK趨勢檢驗對徑流序列進行判斷。線性傾向估計的原理是利用一元線性回歸方程擬合時間序列和觀測值，并以回歸系數作為判斷序列長期趨勢的方法[18]。其具體求解公式為

y=ax+b

(2)

(3)

式中，y為徑流；x為年份；a為回歸系數(線性傾向估計值)。

MK趨勢檢驗屬于非參數檢驗方法，因其不需要樣本遵循某一特定的分布，并且不易受到異常值的干擾，常用于氣象水文要素時間序列的趨勢變化檢驗[19]。階段和突變點的檢驗利用累積距平和MK突變檢驗方法進行[18]，徑流序列分解各模態的周期利用功率譜分析方法進行識別[20]。

2 贛江徑流變化特征

2.1 年際變化與趨勢

根據贛江外洲站1950年～2015年，共計66 a的徑流觀測資料，繪制出徑流變化序列(見圖1)。

圖1 1950年～2015年贛江徑流變化序列

由圖1可以發現，贛江年徑流變化波動范圍較大，最小值出現于1963年，徑流量為236.07億m3；最大值出現于1973年，為1147.87億m3；其多年平均值為682.46億m3。線性傾向估計的結果表明，贛江徑流變化具有上升趨勢，其上升速率為6.5×108m3/10 a。MK趨勢檢驗的結果表明，其Z統計量為0.653。由于其Z值<1.64，因此未能通過相關顯著性檢驗，也就是說從長期來看贛江徑流變化的上升趨勢并不顯著。

累積距平分析的結果顯示(見圖2)，近66年來贛江徑流變化可以分為4個階段：第一階段為1950年～1969年，這一時間段內流量呈波動下降趨勢。第二階段為1970年～1990年，這一時間段流量變化呈振動變化，變化較為劇烈，但總體屬于流量較少的時期。第三階段為1991年～2002年，這一時間段徑流呈明顯增加的趨勢。第四階段為2003年～2015年，這一時間段流量較高，但呈下降趨勢。進一步的MK突變檢驗表明(見圖3)，贛江徑流序列在1952年～1954年、1989年～1990年、2003年～2014年3個時間段發生了突變，由于序列的兩端可能會受到邊界效應的影響，其檢測可能存在一定的誤差，一般認為處于序列中間的突變點比較符合實際情況[20]。由圖3可以觀測到：1989年/1990年UF統計量和UB統計量有交點，并且處于±1.96范圍之內，這說明這一突變點通過了α=0.05的顯著性檢驗。將整個贛江徑流序列從這一突變點分開，分別統計前后兩個時間段的徑流平均值可以發現，1950年～1989年的平均值為663.47億m3，1990年～2015年的平均值為711.69億m3。

圖2 1950年～2015年徑流序列累積距平

圖3 1950年～2015年徑流序列MK突變檢驗

2.2 多尺度變化特征

2.2.1 模態特征

利用CEEMD方法對1950年～2015年贛江徑流變化序列進行分解(見圖4)，共得到5個IMF分量和一個趨勢分量(RES)。其中，每個IMF分量對應于一個窄波段的徑流變化波動，反映特定時間尺度下徑流量的周期變化和振幅隨時間的演變過程[21]。長期趨勢項代表的是序列中存在的緩慢變化趨勢，其變化周期大于原始序列的總長度。Huang等[22]的研究指出，如果趨勢項為單調上升或單調下降函數則說明原始信號具有明顯的非平穩性。由圖4可以看出RES項是單調上升，因此也證實了贛江徑流變化序列屬于非平穩性的復雜變化過程。

圖4 1950年～2015年徑流序列的CEEMD分解結果

由圖3可觀察到，序列的各IMF分量變化規律符合自然信號的非線性變化的特征，都具有相對穩定的變化周期。不過各IMF分量對原始序列的影響程度有差異，IMF1的規律性較差，仍然具有一定的無序變化過程。但是IMF2～IMF5已經具有比較明顯的周期變化規律。相關研究已經證明各IMF分量具有獨立的代表性，不會在同一時間出現周期相同的振蕩。而且各IMF分量一般具有物理意義，揭示了原始信號中固有的不同尺度的波動過程，反映了氣象水文復雜系統內部動力過程和外部強迫共同作用的非線性特征[22]。經過對各IMF和RES項的平均周期和振幅分析，得到分解結果的主要統計值(見表1)。IMF1～IMF5的平均周期在3.1～66 a范圍內，其最大振幅隨著分解的遞進而下降，各分量的方差解釋量在63.6%至0.29%之間，其中前4項的總方差解釋量為99.5%。各分量與原始序列的相關分析表明，前4項的相關系數較高，在0.79至0.24之間，并通過了α=0.01顯著性檢驗。

表1 IMF分量和趨勢分量的主要統計值

2.2.2 周期特征

為了進一步揭示出不同時間尺度下贛江徑流變化的周期特征，同時結合上一節分析結果中所得出的CEEMD分解前4項方差之和占到序列總方差的99.5%這一特征，利用功率譜分析方法對IMF1～IMF4進行周期分析，以確定其所包含的全部周期信息。結果顯示(見圖5)，IMF1具有3.6 、4.2 、4.4 a周期；IMF2具有3.8、4.4、4.7、7.3、8.0 a周期；IMF3具有8.9、10、11.4、13.3、16、20 a周期；IMF4具有16、20、26.7、40 a周期。

圖5 各IMF分量的功率譜分析(虛線表示95%置信區間)

3 討 論

CEEMD分解的結果表明，贛江徑流變化具有明顯的非線性特征，各IMF分量代表了不同時間尺度下的徑流變化規律。根據表1中的統計結果可以發現，IMF1的方差解釋量和相關系數都是最大的，因此可以判斷出在年際尺度上序列具有顯著的準3 a變化周期。同時，IMF1的功率譜分析的結果也印證了其具有3～4 a的年際變化周期。黃榮輝等[23]利用現代降水資料的分析表明，長江流域降水具有3 a左右的周期變化特征。對長江流域梅雨降水量的分析也表明，準3 a周期是其年際變化的主要表現形式[24]。由圖4可以觀察到，除1976年～1990年這一段時間IMF1的振幅較小外，其余時間段振蕩信號十分明顯，而且自2000年以來有振幅增大的趨勢。功率譜分析表明IMF2具有準4～8 a的變化周期，結合表1中的5.7 a平均周期，可以判斷其具有準6 a 的周期變化規律。IMF2方差解釋量為16.4%，代表了贛江徑流變化的次顯著變化特征，其振幅在1950年～1960年較大，其后有逐漸減小的趨勢。王鐘睿和錢永甫的研究表明[25]，4～8 a周期與厄爾尼諾有關；胡光偉等[26]對湘江徑流的周期分析也表明，準6 a周期表現明顯。

在年代際尺度上，IMF3和IMF4的平均周期為準14 a和準26 a。其功率譜分析也同樣表明，準14 a和準26 a的譜值最大，信號最為顯著?？紤]到IMF3和IMF4的方差解釋量只占到總量的10%左右，因此其體現的振蕩過程所包含的實際物理意義可能相對較小。相關研究表明[27]，太陽活動是在8～15 a之間變化，而且26 a周期也與太陽活動的長期波動相關。因此贛江徑流變化的準14 a和準26 a周期很可能是太陽活動這一外部強迫的體現。另一方面，IMF5和RES趨勢項所體現的是序列較長時間尺度下的變化趨勢，結合圖4可知其未來發展有增加的趨勢。但是考慮到MK趨勢檢驗的結果，這一增加的趨勢可能并不顯著。

4 結 論

CEEMD分解是一種適用于非線性復雜序列的信號分解方法，可以提取原始信號背后的不同尺度的變化過程。本文利用CEEMD方法，結合趨勢分析、突變檢驗和功率譜分析等手段，對1950年～2015年贛江外洲站徑流序列進行多尺度變化特征分析，得到以下主要結論：

(1)線性傾向估計的結果表明，贛江徑流變化上升速率為6.5×108m3/10 a。MK趨勢檢驗的結果表明，其Z統計量為0.653，上升趨勢并不顯著。

(2)累積距平分析的結果顯示，近66年來贛江徑流變化可以分為4個階段。其中下降階段為1950年～1969年、2003年～2015年；上升階段為1991年～2002年；波動階段為1970年～1990年。MK突變檢驗表明1989年/1990年為突變點，并通過了α=0.05的顯著性檢驗。

(3)CEEMD的分解結果表明，贛江徑流變化可以被分解為5個IMF分量和一個趨勢分量(RES)。IMF1～IMF5的平均周期在3.1～66 a范圍內，其最大振幅隨著分解的遞進而下降，前4項的總方差解釋量為99.5%。各分量與原始序列的相關分析表明，前4項的相關系數較高，在0.79至0.24之間，并通過了α=0.01顯著性檢驗。

(4)結合CEEMD和功率譜分析結果表明，在年際尺度上，贛江徑流變化具有準3 a和準6 a周期；在年代際尺度上，贛江徑流變化具有準14 a和準26 a 周期。徑流的準3 a周期可能與長江流域降水的準3 a周期相關，準6 a、準14 a和準26 a周期則可能與厄爾尼諾和太陽活動相關。