基于ESMD的西江流域降水時空變化特征分析

李繼清，陳思雨

（華北電力大學水利與水電工程學院，北京102206）

0 引 言

氣候與人類社會生活密切相關，隨著全球變暖問題日益嚴重［1］，進一步加劇了降水這一氣候基本因子的時空分布差異，進而對水資源、生態環境系統以及人類的生產生活產生影響。西江流域擁有豐富的水資源，降水是其徑流的補給來源，決定徑流分布情況。研究西江流域降水中長期的時空變化特性，能夠掌握該流域水文資源的周期、趨勢和突變規律，為該地區的水資源利用和管理提供科學依據。

氣候由多個時空尺度構成了多層次結構，是一個具有非平穩性的非線性復雜系統。氣候變化特性研究領域已有較為系統的分析方法，所采用的數學分析方法也在逐步改進和完善。早期研究中多采用累積距平、等值線等方法［2-5］，后期逐漸引入數學分析法，國內外在研究各地區降水變化時廣泛采用Mann-Kendall 檢驗法［6，7］，并結合小波分析等多種方法進行對比［8-12］。此外，有序聚類、全球氣候模型等方法也可應用于氣候分析［13-17］。但傳統方法受其理論基礎限制，在實際應用中存在一定的局限性，極點對稱模態分解（Extreme-point Symmetric Mode Decomposition，ESMD）法作為近年來發展起來的局部自適應時間序列分析技術［18］，善于從非平穩、非線性的觀測序列中尋找變化趨勢并對其做出異常診斷。王金良和李宗軍［19］以及房賢水［20］分別從ESMD 方法統計特征、工作原理等方面對其優勢進行了論證。近年來，該方法也逐漸應用到氣候變化的實例研究［21-23］，進一步驗證了其優勢和可靠性。本文采用ESMD 法，對西江流域的降水觀測資料進行周期、突變及趨勢分析，并與線性趨勢分析、Mann-Kendall 檢驗法和小波分析法的分析結果進行對比驗證，進而詳細分析西江流域的降水時空變化特征，為流域水資源的開發利用提供決策依據。

1 研究方法

極點對稱模態分解法（ESMD）能夠通過模態分解和時-頻分析，全面分析時間序列的周期、突變和趨勢特征，在水文時間序列分析領域具有較好的優勢。分別利用線性趨勢分析法判別數據序列的整體線性變化趨勢，Mann-Kendall 檢驗法判別數據數列的變化趨勢和突變時間點，小波分析法提取數據序列的周期特征，與ESMD方法的計算結果進行對比分析。

1.1 極點對稱模態分解法（ESMD）

極點對稱模態分解（ESMD）方法在希爾伯特-黃變換法之上發展而來，在保留了經驗模態分解法優點的同時也解決了其中存在的“模態混疊問題”，該方法在氣候分析中的優勢為：擅長尋找變化趨勢，能夠分離出數據序列的年際變化趨勢和整體趨勢；擅長異常診斷，能夠觀察數據序列各模態中的異常時頻段；擅長時-頻分析，采用直接插值法能夠更為直觀地分析各時間尺度上的頻率變化［24］。

第一部分模態分解具體計算步驟如下：

（1）找出序列X的全部極值點，包括極大值點和極小值點，并將它們依次記作Ei（i=1，2，3，…，n）；

（2）用線段連接所有相鄰的極點并將線段中點依次記為Fi（1≤i≤n-1）；并通過一定方式補充左、右兩端邊界中點F0和Fn；

（3）利用所獲取的n+1 個中點構建p條內插曲線L1，…，Lp（p≥1），并計算其均值曲線L*；

（4）對X-L*序列重復上述步驟，直到|L*|≤ε 或者篩選次數達到預設的最大值K，得到第一個經驗模M1；

（5）對剩余序列X-M1重復上述4 個步驟，直到剩余序列R只剩余一定數量的極點，便可分別得到經驗模M2，M3，…；

（6）讓最大篩選次數K在整數區間［Kmin，Kmax］內變化并重復上述幾個步驟。計算方差比率，并畫出隨K的變化圖，找出方差比率最小值對應的最大篩選次數K0，此時分解得到的趨勢余項R為數據的最佳擬合曲線。以K0作為限制條件再次重復上述步驟，得到最佳分解結果。最終利用ESMD 方法將序列X分解為一系列經驗模態和一個趨勢余項。

第二部分時-頻分析認為應當尊重離散數據的離散性特征，提出了針對數據的“直接插值法”，不但可以直觀地體現各模態的振幅與頻率的時變性，還可明確地獲知總能量變化。直接插值法基本思路如下：

（1）尋找極值點，計算兩個相鄰極大值點和相鄰極小值點之間的時間差；

（2）將這些時間段視為局部周期賦給其中點，畫出時間-周期對應點圖；

（3）將這些局部周期值取倒數得到局部頻率，再做3次樣條插值得到光滑的時間-頻率變化曲線［25］。

圖1 ESMD方法流程圖Fig.1 The step of ESMD

1.2 對比驗證方法

1.2.1 線性趨勢法

線性趨勢法能夠對數據序列的線性變化趨勢進行分析。用yi表示樣本量為n的氣候變量，ti為對應時刻，建立一元線性回歸方程：

式中：a為回歸常數；b為回歸系數即氣候變量的傾向趨勢［26］。

1.2.2 Mann-Kendall檢驗法

Mann-Kendall 檢驗法是一種非參數統計檢驗方法，廣泛應用于水文序列的趨勢、突變特征分析。利用M-K 統計值Z進行趨勢統計的顯著性檢驗，先假設該序列無趨勢，通過兩尾檢驗，在給定顯著性水平下，查得臨界值Z1-α/2，當|Z|＜Z1-α/2時，接受原假設；若|Z|＞Z1-α/2，則拒絕原假設。進行突變分析需要計算M-K統計值UFk和UBk：

通過分析UFk和UBk交點位置即可明確突變發生時間［27］。

1.2.3 小波分析法

小波分析的基本思想是用一簇小波函數系來表示或逼近某一信號或函數。在小波變換系數Wf（c，d）中，將小波變換系數的平方值在d域上積分可得到小波方差：

式中：c為伸縮尺度；d為平移參數。

小波方差隨尺度c的變化過程為小波方差圖，能夠反應信號波動的能量隨尺度c的分布，進而獲得數據序列的周期變化規律［28］。

2 研究結果與分析

采用ESMD 法對于西江流域的降水資料進行模態分解和時-頻分析，結合線性趨勢分析、Mann-Kendall 檢驗法和小波分析法，對各降水系列的周期、突變和趨勢特征進行分析討論。

2.1 研究區概況及數據資料

西江是珠江流域的主干河流，發源于云南省，流經地區地勢西北高東南低，全長共2 075 km，主要河流由南盤江、紅水河等河段組成。西江流域屬于亞熱帶氣候區，春夏季降雨充沛，秋冬較為干旱，水資源量變化主要受到降水季節分配的影響，徑流年內分配極不均勻，汛期徑流量比例占年徑流量的80%以上。降水量由東到西呈現遞減趨勢［29］。

分析數據來自于西江流域1954-2008年51 個氣象站的日降水資料，空間分布如圖2，基本能夠均勻覆蓋西江整個流域。對日降水資料進行整理計算，得到西江流域年降水、年降水極大值、年降水極小值、汛期降水和非汛期降水數據系列。其中年降水極值為全年日降水極大或極小值。根據地理情況將西江流域劃分為500 m 以上的上游高原區以及500 m 以下的中下游平原盆地區。平原區年降水Cv值為0.121，高原區年降水Cv值為0.089，圖3為兩地區年降水量對比。高原區年降水量約為1 100 mm，平原區年降水量約為1 500 mm，兩地區降水趨勢基本相同。

圖2 西江流域氣象站分布圖Fig.2 Distribution of meteorological stations in the Xijiang River Basin

圖3 西江流域平原區和高原區年降水概況圖Fig.3 Annual precipitation in the plain and plateau of the Xijiang River Basin

2.2 周期性

ESMD 法具有自適應性和基于信號的局部變化特性，能夠更好地提取非平穩、非線性序列的變化規律。通過ESMD 法對西江流域各降水序列進行模態分解，可得到3～5 個模態分量以及一個趨勢項余量R，玉溪和連縣年降水分解結果圖4，各模態分量能夠反應原數據序列中固有的振蕩特征，即對應為降水序列的周期規律。表1為各降水序列的周期結果。西江流域各站之間各降水序列的固定周期成分均出現接近或相同的情況，具有統一的空間周期特征。

圖4 年降水模態分解圖Fig.4 Exploded modal decomposition maps of annual precipitation

表1 西江流域內各站各降水指標周期統計表 aTab.1 Periodic statistics of precipitation indicators at stations in the Xijiang River Basin

對西江流域各站各降水指標周期成分出現次數進行統計整理，結果如圖5。其中橙色代表年降水，紫色代表年降水極大值，綠色代表年降水極小值，黃色代表汛期，藍色代表非汛期。年降水18 a 的周期規律在站中出現頻率最高，其次為28 a 的周期特征。年極大值和年極小值的周期規律出現頻次最高的均為28 a，其次為2 a、14 a 等。汛期降水以存在3 a 和28 a 的周期特征最為普遍。非汛期降水序列中84%的氣象站存在2 a 的周期變化特征，其次是28 a 和18 a。各降水指標之間的周期特征具有較高一致性，主要存在2、14、18和28 a的周期變化規律。

圖5 各降水指標周期成分出現次數統計圖Fig.5 Statistics of the number of occurrences of the periodic components of various precipitation index

采用小波分析法與ESMD 方法周期分析結果進行對比，表2為貴陽、羅甸、通道和威寧年降水的周期對比結果。兩方法下存在部分一致或相近的周期結果，周期越短結果一致性越好，周期越長越容易出現差異。究其原因，小波分析是以傅里葉變換為理論依據，在小波基選取等方面存在一定局限性，限制了該方法的實際運用效果。而ESMD方法具有較強的靈活性和自適應性，可通過數據自身特點進行分解，探索事物的內在規律。

表2 ESMD法和小波分析法下年降水周期對比結果表 aTab.2 Comparison of sudden changes in precipitation in the next year by ESMD method and wavelet analysis

2.3 突變性

根據ESMD 法得到的數據序列模態分量，采用直接插值法得到其對應的頻率-振幅時變圖，如圖6為玉溪和連縣站年降水時-頻分析結果圖，其中F代表頻率，A代表振幅，根據兩者對應關系，可知圖像在出現低頻、大幅度振蕩和高頻、小幅度振蕩時均代表數據序列發生了降水突變。對比圖中F和A的變化情況，可知玉溪年降水在1983年和1993年發生，連縣年降水在1971年發生突變。

圖6 年降水頻率-振幅時變圖Fig.6 Frequency-amplitude time-varying data of annual precipitation

通過觀察頻率-振幅時變圖得到各降水序列突變年份，結果整理統計如表3。空間上，各降水指標在站間的突變結果均不完全一致，在相鄰站間突變結果一致性稍高，圖7梧州、廣寧、高要3 個相鄰氣象站的年降水均在觀測期內發生三次突變現象，且突變時間較為接近。但降水突變在西江流域整體空間范圍內未顯示出明顯的空間分布規律和特征。

圖7 梧州、廣寧和高要站的年降水突變結果圖Fig.7 Abrupt changes in annual precipitation at Wuzhou，Guangning and Gaoyao stations

表3 西江流域內各站各降水指標突變年份統計表Tab.3 Statistics of abrupt change of precipitation index of each station in Xijiang River Basin

通過整理各降水指標的突變時間點分布年代，分析突變的時間分布特征，結果如圖8。各降水指標在各年代內均有突變現象發生，其中年降水和汛期降水的突變結果相似，均在20世紀60年代出現突變現象最為頻繁，其次為90年代；非汛期降水在80年代出現突變現象的次數最多，60年代的突變次數與汛期降水持平。年極大值的突變現象頻發于60年代和90年代，年極小值則在80年代出現突變情況最多。各降水指標突變現象多集中在60年代、80年代和90年代，在50年代和00年代較少發生突變。

圖8 各降水指標突變時間出現年代統計圖Fig.8 The chronological statistics of the abrupt time of each precipitation index

將ESMD 法與Mann-Kendall 檢驗法結果進行對比分析，圖9為玉溪和連縣的年降水M-K 結果圖，其中玉溪在1956年、1977年、1993年、2001年和2008年發生降水突變，連縣站的突變年份為1960年、1963年、1970年和1986年。

圖9 年降水M-K圖Fig.9 Results of annual precipitation based on M-K method

表4為蒙自、獨山、欽州和百色站年降水M-K法和ESMD法的突變結果對比。雖然兩方法算得突變結果不完全一致，但均存在相近或相同的突變結果，可知ESMD 方法進行西江流域各降水序列突變分析具有一定的可靠性。兩方法相比，M-K 法僅能對數據序列中的突變點進行判別，而ESMD 方法能夠在通過各模態頻率與振幅的對應關系判斷數據序列突變點的同時了解其局部變化情況，具有明顯優勢。

表4 ESMD法和Mann-Kendall檢驗法下部分站年降水突變年份對比結果表Tab.4 Results of abrupt change in precipitation of some stations under the ESMD method and the M-K test

2.4 趨勢性

線性趨勢分析能夠判別西江流域各降水指標的多年線性變化特征，統計結果如表5。年降水、年極大值和汛期降水線性趨勢分布情況相似，正、負增長站數基本持平，3 個指標的正增長率最大均為防城站，負增長率最大均為南寧站，其中年極大值的線性變化速率略小，多不超過1 mm/a，年降水和汛期降水線性變化速率多不超過4 mm/a。年極小值和非汛期降水的線性趨勢多為正增長，但線性變化趨勢均不顯著，年極小值正增長率最大為連縣站，負增長率最大為融安站，非汛期降水的正增長率最大為靖西站，負增長率最大為高要站。

表5 西江流域各站各降水指標線性變化速率結果表Tab.5 Results of the linear change rate of precipitation indexes at various stations in the Xijiang River Basin

ESMD 方法模態分解得到的趨勢項R可描述降水序列的具體變化過程，進而分析西江流域各降水指標趨勢特征。圖10為年降水指標的變化趨勢類型，主要變化趨勢為上升-下降-上升型和先降后升型，多分布于中下游地區，除此以外上游高原區域存在少量的先升后降型和下降-上升-下降型。年極大值指標的變化趨勢以上升-下降-上升、先升后降以及先降后升型為主，其中前兩種變化類型均多分布于流域中下游，此外上游地區還存在少量持續上升和下降-上升-下降的變化特征。年極小值指標中，持續上升型、先降后升型、先升后降型等多種變化類型數量較為均勻，在流域內交錯分布。汛期降水與年降水分布情況相似，出現了少量站點呈現持續下降的變化特征。非汛期降水指標中近半數站呈現先升后降的變化特點，多分布于流域中下游。由此可見，西江流域中下游地區各降水指標變化趨勢類型較為統一，而位于流域西部的上游地區由于山區地形因素影響，各降水指標的趨勢變化存在較大差異。

圖10 ESMD方法下年降水的各類變化趨勢圖Fig.10 Various trends of annual precipitation based on the ESMD method

采用M-K 趨勢檢驗法對ESMD 法的趨勢結果進行對比分析，M-K 趨勢檢驗法給定顯著水平α=0.05，臨界值為±1.96，兩方法計算得到的桂平站年降水趨勢結果對比如圖11。在M-K法分析中，桂平站年降水前期呈現較為平穩的小幅波動，在90年代后呈現增長趨勢，整體降水序列M-K 統計值Z為1.19，表現為不顯著上升趨勢，與ESMD 法的趨勢結果具有較高的一致性，證明采用ESMD方法分析數據序列的趨勢性具有可靠性。

圖11 桂平站年降水ESMD法與M-K法趨勢結果對比圖Fig.11 Comparison of trend results of ESMD and M-K methods of annual precipitation at Guiping station

3 結 論

氣候變化與水文資源有著密切關聯，通過研究降水變化特性能為流域水資源管理提供有效手段。在氣候分析領域，ESMD 方法的優勢在于能夠更為全面有效地識別數據序列的趨勢、周期和突變特征。采用ESMD 方法結合小波分析等對比方法分析西江流域的各降水指標變化特征，得到如下結論：①西江流域各站的降水指標周期特征具有統一性，主要存在2、14、18、28 a等周期規律；②每個站的各降水序列均存在1～5個突變點不等，各降水指標突變多集中在20世紀60年代、80年代和90年代；③年降水、年極大值和汛期降水的正、負增長線性趨勢站數量基本持平，年極小值和非汛期降水的線性趨勢多呈現出上升現象；④年降水和非汛期降水在中下游區域的ESMD 法趨勢變化具有較高的一致性，其余各降水指標的變化趨勢多呈現散亂分布。

在全球變暖環境下，西江流域降水指標年內、年際周期變化規律相近，年極小值及非汛期降水量統一呈現微弱增長趨勢，中下游地區降水指標趨勢特征相對一致，降水偏多和偏少時期交替出現，但由于地理位置等因素影響，各站點仍存在其獨有的突變和趨勢特征，需進一步探尋各水文要素的變化特性及各要素之間的相關聯系和影響機制。