李繼清,陳思雨
(華北電力大學水利與水電工程學院,北京102206)
氣候與人類社會生活密切相關,隨著全球變暖問題日益嚴重[1],進一步加劇了降水這一氣候基本因子的時空分布差異,進而對水資源、生態(tài)環(huán)境系統(tǒng)以及人類的生產(chǎn)生活產(chǎn)生影響。西江流域擁有豐富的水資源,降水是其徑流的補給來源,決定徑流分布情況。研究西江流域降水中長期的時空變化特性,能夠掌握該流域水文資源的周期、趨勢和突變規(guī)律,為該地區(qū)的水資源利用和管理提供科學依據(jù)。
氣候由多個時空尺度構(gòu)成了多層次結(jié)構(gòu),是一個具有非平穩(wěn)性的非線性復雜系統(tǒng)。氣候變化特性研究領域已有較為系統(tǒng)的分析方法,所采用的數(shù)學分析方法也在逐步改進和完善。早期研究中多采用累積距平、等值線等方法[2-5],后期逐漸引入數(shù)學分析法,國內(nèi)外在研究各地區(qū)降水變化時廣泛采用Mann-Kendall 檢驗法[6,7],并結(jié)合小波分析等多種方法進行對比[8-12]。此外,有序聚類、全球氣候模型等方法也可應用于氣候分析[13-17]。但傳統(tǒng)方法受其理論基礎限制,在實際應用中存在一定的局限性,極點對稱模態(tài)分解(Extreme-point Symmetric Mode Decomposition,ESMD)法作為近年來發(fā)展起來的局部自適應時間序列分析技術(shù)[18],善于從非平穩(wěn)、非線性的觀測序列中尋找變化趨勢并對其做出異常診斷。王金良和李宗軍[19]以及房賢水[20]分別從ESMD 方法統(tǒng)計特征、工作原理等方面對其優(yōu)勢進行了論證。近年來,該方法也逐漸應用到氣候變化的實例研究[21-23],進一步驗證了其優(yōu)勢和可靠性。本文采用ESMD 法,對西江流域的降水觀測資料進行周期、突變及趨勢分析,并與線性趨勢分析、Mann-Kendall 檢驗法和小波分析法的分析結(jié)果進行對比驗證,進而詳細分析西江流域的降水時空變化特征,為流域水資源的開發(fā)利用提供決策依據(jù)。
極點對稱模態(tài)分解法(ESMD)能夠通過模態(tài)分解和時-頻分析,全面分析時間序列的周期、突變和趨勢特征,在水文時間序列分析領域具有較好的優(yōu)勢。分別利用線性趨勢分析法判別數(shù)據(jù)序列的整體線性變化趨勢,Mann-Kendall 檢驗法判別數(shù)據(jù)數(shù)列的變化趨勢和突變時間點,小波分析法提取數(shù)據(jù)序列的周期特征,與ESMD方法的計算結(jié)果進行對比分析。
極點對稱模態(tài)分解(ESMD)方法在希爾伯特-黃變換法之上發(fā)展而來,在保留了經(jīng)驗模態(tài)分解法優(yōu)點的同時也解決了其中存在的“模態(tài)混疊問題”,該方法在氣候分析中的優(yōu)勢為:擅長尋找變化趨勢,能夠分離出數(shù)據(jù)序列的年際變化趨勢和整體趨勢;擅長異常診斷,能夠觀察數(shù)據(jù)序列各模態(tài)中的異常時頻段;擅長時-頻分析,采用直接插值法能夠更為直觀地分析各時間尺度上的頻率變化[24]。
第一部分模態(tài)分解具體計算步驟如下:
(1)找出序列X的全部極值點,包括極大值點和極小值點,并將它們依次記作Ei(i=1,2,3,…,n);
(2)用線段連接所有相鄰的極點并將線段中點依次記為Fi(1≤i≤n-1);并通過一定方式補充左、右兩端邊界中點F0和Fn;
(3)利用所獲取的n+1 個中點構(gòu)建p條內(nèi)插曲線L1,…,Lp(p≥1),并計算其均值曲線L*;
(4)對X-L*序列重復上述步驟,直到|L*|≤ε 或者篩選次數(shù)達到預設的最大值K,得到第一個經(jīng)驗模M1;
(5)對剩余序列X-M1重復上述4 個步驟,直到剩余序列R只剩余一定數(shù)量的極點,便可分別得到經(jīng)驗模M2,M3,…;
(6)讓最大篩選次數(shù)K在整數(shù)區(qū)間[Kmin,Kmax]內(nèi)變化并重復上述幾個步驟。計算方差比率,并畫出隨K的變化圖,找出方差比率最小值對應的最大篩選次數(shù)K0,此時分解得到的趨勢余項R為數(shù)據(jù)的最佳擬合曲線。以K0作為限制條件再次重復上述步驟,得到最佳分解結(jié)果。最終利用ESMD 方法將序列X分解為一系列經(jīng)驗模態(tài)和一個趨勢余項。
第二部分時-頻分析認為應當尊重離散數(shù)據(jù)的離散性特征,提出了針對數(shù)據(jù)的“直接插值法”,不但可以直觀地體現(xiàn)各模態(tài)的振幅與頻率的時變性,還可明確地獲知總能量變化。直接插值法基本思路如下:
(1)尋找極值點,計算兩個相鄰極大值點和相鄰極小值點之間的時間差;
(2)將這些時間段視為局部周期賦給其中點,畫出時間-周期對應點圖;
(3)將這些局部周期值取倒數(shù)得到局部頻率,再做3次樣條插值得到光滑的時間-頻率變化曲線[25]。
圖1 ESMD方法流程圖Fig.1 The step of ESMD
1.2.1 線性趨勢法
線性趨勢法能夠?qū)?shù)據(jù)序列的線性變化趨勢進行分析。用yi表示樣本量為n的氣候變量,ti為對應時刻,建立一元線性回歸方程:
式中:a為回歸常數(shù);b為回歸系數(shù)即氣候變量的傾向趨勢[26]。
1.2.2 Mann-Kendall檢驗法
Mann-Kendall 檢驗法是一種非參數(shù)統(tǒng)計檢驗方法,廣泛應用于水文序列的趨勢、突變特征分析。利用M-K 統(tǒng)計值Z進行趨勢統(tǒng)計的顯著性檢驗,先假設該序列無趨勢,通過兩尾檢驗,在給定顯著性水平下,查得臨界值Z1-α/2,當|Z|<Z1-α/2時,接受原假設;若|Z|>Z1-α/2,則拒絕原假設。進行突變分析需要計算M-K統(tǒng)計值UFk和UBk:
通過分析UFk和UBk交點位置即可明確突變發(fā)生時間[27]。
1.2.3 小波分析法
小波分析的基本思想是用一簇小波函數(shù)系來表示或逼近某一信號或函數(shù)。在小波變換系數(shù)Wf(c,d)中,將小波變換系數(shù)的平方值在d域上積分可得到小波方差:
式中:c為伸縮尺度;d為平移參數(shù)。
小波方差隨尺度c的變化過程為小波方差圖,能夠反應信號波動的能量隨尺度c的分布,進而獲得數(shù)據(jù)序列的周期變化規(guī)律[28]。
采用ESMD 法對于西江流域的降水資料進行模態(tài)分解和時-頻分析,結(jié)合線性趨勢分析、Mann-Kendall 檢驗法和小波分析法,對各降水系列的周期、突變和趨勢特征進行分析討論。
西江是珠江流域的主干河流,發(fā)源于云南省,流經(jīng)地區(qū)地勢西北高東南低,全長共2 075 km,主要河流由南盤江、紅水河等河段組成。西江流域?qū)儆趤啛釒夂騾^(qū),春夏季降雨充沛,秋冬較為干旱,水資源量變化主要受到降水季節(jié)分配的影響,徑流年內(nèi)分配極不均勻,汛期徑流量比例占年徑流量的80%以上。降水量由東到西呈現(xiàn)遞減趨勢[29]。
分析數(shù)據(jù)來自于西江流域1954-2008年51 個氣象站的日降水資料,空間分布如圖2,基本能夠均勻覆蓋西江整個流域。對日降水資料進行整理計算,得到西江流域年降水、年降水極大值、年降水極小值、汛期降水和非汛期降水數(shù)據(jù)系列。其中年降水極值為全年日降水極大或極小值。根據(jù)地理情況將西江流域劃分為500 m 以上的上游高原區(qū)以及500 m 以下的中下游平原盆地區(qū)。平原區(qū)年降水Cv值為0.121,高原區(qū)年降水Cv值為0.089,圖3為兩地區(qū)年降水量對比。高原區(qū)年降水量約為1 100 mm,平原區(qū)年降水量約為1 500 mm,兩地區(qū)降水趨勢基本相同。
圖2 西江流域氣象站分布圖Fig.2 Distribution of meteorological stations in the Xijiang River Basin
圖3 西江流域平原區(qū)和高原區(qū)年降水概況圖Fig.3 Annual precipitation in the plain and plateau of the Xijiang River Basin
ESMD 法具有自適應性和基于信號的局部變化特性,能夠更好地提取非平穩(wěn)、非線性序列的變化規(guī)律。通過ESMD 法對西江流域各降水序列進行模態(tài)分解,可得到3~5 個模態(tài)分量以及一個趨勢項余量R,玉溪和連縣年降水分解結(jié)果圖4,各模態(tài)分量能夠反應原數(shù)據(jù)序列中固有的振蕩特征,即對應為降水序列的周期規(guī)律。表1為各降水序列的周期結(jié)果。西江流域各站之間各降水序列的固定周期成分均出現(xiàn)接近或相同的情況,具有統(tǒng)一的空間周期特征。
圖4 年降水模態(tài)分解圖Fig.4 Exploded modal decomposition maps of annual precipitation
表1 西江流域內(nèi)各站各降水指標周期統(tǒng)計表 aTab.1 Periodic statistics of precipitation indicators at stations in the Xijiang River Basin
對西江流域各站各降水指標周期成分出現(xiàn)次數(shù)進行統(tǒng)計整理,結(jié)果如圖5。其中橙色代表年降水,紫色代表年降水極大值,綠色代表年降水極小值,黃色代表汛期,藍色代表非汛期。年降水18 a 的周期規(guī)律在站中出現(xiàn)頻率最高,其次為28 a 的周期特征。年極大值和年極小值的周期規(guī)律出現(xiàn)頻次最高的均為28 a,其次為2 a、14 a 等。汛期降水以存在3 a 和28 a 的周期特征最為普遍。非汛期降水序列中84%的氣象站存在2 a 的周期變化特征,其次是28 a 和18 a。各降水指標之間的周期特征具有較高一致性,主要存在2、14、18和28 a的周期變化規(guī)律。
圖5 各降水指標周期成分出現(xiàn)次數(shù)統(tǒng)計圖Fig.5 Statistics of the number of occurrences of the periodic components of various precipitation index
采用小波分析法與ESMD 方法周期分析結(jié)果進行對比,表2為貴陽、羅甸、通道和威寧年降水的周期對比結(jié)果。兩方法下存在部分一致或相近的周期結(jié)果,周期越短結(jié)果一致性越好,周期越長越容易出現(xiàn)差異。究其原因,小波分析是以傅里葉變換為理論依據(jù),在小波基選取等方面存在一定局限性,限制了該方法的實際運用效果。而ESMD方法具有較強的靈活性和自適應性,可通過數(shù)據(jù)自身特點進行分解,探索事物的內(nèi)在規(guī)律。
表2 ESMD法和小波分析法下年降水周期對比結(jié)果表 aTab.2 Comparison of sudden changes in precipitation in the next year by ESMD method and wavelet analysis
根據(jù)ESMD 法得到的數(shù)據(jù)序列模態(tài)分量,采用直接插值法得到其對應的頻率-振幅時變圖,如圖6為玉溪和連縣站年降水時-頻分析結(jié)果圖,其中F代表頻率,A代表振幅,根據(jù)兩者對應關系,可知圖像在出現(xiàn)低頻、大幅度振蕩和高頻、小幅度振蕩時均代表數(shù)據(jù)序列發(fā)生了降水突變。對比圖中F和A的變化情況,可知玉溪年降水在1983年和1993年發(fā)生,連縣年降水在1971年發(fā)生突變。
圖6 年降水頻率-振幅時變圖Fig.6 Frequency-amplitude time-varying data of annual precipitation
通過觀察頻率-振幅時變圖得到各降水序列突變年份,結(jié)果整理統(tǒng)計如表3??臻g上,各降水指標在站間的突變結(jié)果均不完全一致,在相鄰站間突變結(jié)果一致性稍高,圖7梧州、廣寧、高要3 個相鄰氣象站的年降水均在觀測期內(nèi)發(fā)生三次突變現(xiàn)象,且突變時間較為接近。但降水突變在西江流域整體空間范圍內(nèi)未顯示出明顯的空間分布規(guī)律和特征。
圖7 梧州、廣寧和高要站的年降水突變結(jié)果圖Fig.7 Abrupt changes in annual precipitation at Wuzhou,Guangning and Gaoyao stations
表3 西江流域內(nèi)各站各降水指標突變年份統(tǒng)計表Tab.3 Statistics of abrupt change of precipitation index of each station in Xijiang River Basin
通過整理各降水指標的突變時間點分布年代,分析突變的時間分布特征,結(jié)果如圖8。各降水指標在各年代內(nèi)均有突變現(xiàn)象發(fā)生,其中年降水和汛期降水的突變結(jié)果相似,均在20世紀60年代出現(xiàn)突變現(xiàn)象最為頻繁,其次為90年代;非汛期降水在80年代出現(xiàn)突變現(xiàn)象的次數(shù)最多,60年代的突變次數(shù)與汛期降水持平。年極大值的突變現(xiàn)象頻發(fā)于60年代和90年代,年極小值則在80年代出現(xiàn)突變情況最多。各降水指標突變現(xiàn)象多集中在60年代、80年代和90年代,在50年代和00年代較少發(fā)生突變。
圖8 各降水指標突變時間出現(xiàn)年代統(tǒng)計圖Fig.8 The chronological statistics of the abrupt time of each precipitation index
將ESMD 法與Mann-Kendall 檢驗法結(jié)果進行對比分析,圖9為玉溪和連縣的年降水M-K 結(jié)果圖,其中玉溪在1956年、1977年、1993年、2001年和2008年發(fā)生降水突變,連縣站的突變年份為1960年、1963年、1970年和1986年。
圖9 年降水M-K圖Fig.9 Results of annual precipitation based on M-K method
表4為蒙自、獨山、欽州和百色站年降水M-K法和ESMD法的突變結(jié)果對比。雖然兩方法算得突變結(jié)果不完全一致,但均存在相近或相同的突變結(jié)果,可知ESMD 方法進行西江流域各降水序列突變分析具有一定的可靠性。兩方法相比,M-K 法僅能對數(shù)據(jù)序列中的突變點進行判別,而ESMD 方法能夠在通過各模態(tài)頻率與振幅的對應關系判斷數(shù)據(jù)序列突變點的同時了解其局部變化情況,具有明顯優(yōu)勢。
表4 ESMD法和Mann-Kendall檢驗法下部分站年降水突變年份對比結(jié)果表Tab.4 Results of abrupt change in precipitation of some stations under the ESMD method and the M-K test
線性趨勢分析能夠判別西江流域各降水指標的多年線性變化特征,統(tǒng)計結(jié)果如表5。年降水、年極大值和汛期降水線性趨勢分布情況相似,正、負增長站數(shù)基本持平,3 個指標的正增長率最大均為防城站,負增長率最大均為南寧站,其中年極大值的線性變化速率略小,多不超過1 mm/a,年降水和汛期降水線性變化速率多不超過4 mm/a。年極小值和非汛期降水的線性趨勢多為正增長,但線性變化趨勢均不顯著,年極小值正增長率最大為連縣站,負增長率最大為融安站,非汛期降水的正增長率最大為靖西站,負增長率最大為高要站。
表5 西江流域各站各降水指標線性變化速率結(jié)果表Tab.5 Results of the linear change rate of precipitation indexes at various stations in the Xijiang River Basin
ESMD 方法模態(tài)分解得到的趨勢項R可描述降水序列的具體變化過程,進而分析西江流域各降水指標趨勢特征。圖10為年降水指標的變化趨勢類型,主要變化趨勢為上升-下降-上升型和先降后升型,多分布于中下游地區(qū),除此以外上游高原區(qū)域存在少量的先升后降型和下降-上升-下降型。年極大值指標的變化趨勢以上升-下降-上升、先升后降以及先降后升型為主,其中前兩種變化類型均多分布于流域中下游,此外上游地區(qū)還存在少量持續(xù)上升和下降-上升-下降的變化特征。年極小值指標中,持續(xù)上升型、先降后升型、先升后降型等多種變化類型數(shù)量較為均勻,在流域內(nèi)交錯分布。汛期降水與年降水分布情況相似,出現(xiàn)了少量站點呈現(xiàn)持續(xù)下降的變化特征。非汛期降水指標中近半數(shù)站呈現(xiàn)先升后降的變化特點,多分布于流域中下游。由此可見,西江流域中下游地區(qū)各降水指標變化趨勢類型較為統(tǒng)一,而位于流域西部的上游地區(qū)由于山區(qū)地形因素影響,各降水指標的趨勢變化存在較大差異。
圖10 ESMD方法下年降水的各類變化趨勢圖Fig.10 Various trends of annual precipitation based on the ESMD method
采用M-K 趨勢檢驗法對ESMD 法的趨勢結(jié)果進行對比分析,M-K 趨勢檢驗法給定顯著水平α=0.05,臨界值為±1.96,兩方法計算得到的桂平站年降水趨勢結(jié)果對比如圖11。在M-K法分析中,桂平站年降水前期呈現(xiàn)較為平穩(wěn)的小幅波動,在90年代后呈現(xiàn)增長趨勢,整體降水序列M-K 統(tǒng)計值Z為1.19,表現(xiàn)為不顯著上升趨勢,與ESMD 法的趨勢結(jié)果具有較高的一致性,證明采用ESMD方法分析數(shù)據(jù)序列的趨勢性具有可靠性。
圖11 桂平站年降水ESMD法與M-K法趨勢結(jié)果對比圖Fig.11 Comparison of trend results of ESMD and M-K methods of annual precipitation at Guiping station
氣候變化與水文資源有著密切關聯(lián),通過研究降水變化特性能為流域水資源管理提供有效手段。在氣候分析領域,ESMD 方法的優(yōu)勢在于能夠更為全面有效地識別數(shù)據(jù)序列的趨勢、周期和突變特征。采用ESMD 方法結(jié)合小波分析等對比方法分析西江流域的各降水指標變化特征,得到如下結(jié)論:①西江流域各站的降水指標周期特征具有統(tǒng)一性,主要存在2、14、18、28 a等周期規(guī)律;②每個站的各降水序列均存在1~5個突變點不等,各降水指標突變多集中在20世紀60年代、80年代和90年代;③年降水、年極大值和汛期降水的正、負增長線性趨勢站數(shù)量基本持平,年極小值和非汛期降水的線性趨勢多呈現(xiàn)出上升現(xiàn)象;④年降水和非汛期降水在中下游區(qū)域的ESMD 法趨勢變化具有較高的一致性,其余各降水指標的變化趨勢多呈現(xiàn)散亂分布。
在全球變暖環(huán)境下,西江流域降水指標年內(nèi)、年際周期變化規(guī)律相近,年極小值及非汛期降水量統(tǒng)一呈現(xiàn)微弱增長趨勢,中下游地區(qū)降水指標趨勢特征相對一致,降水偏多和偏少時期交替出現(xiàn),但由于地理位置等因素影響,各站點仍存在其獨有的突變和趨勢特征,需進一步探尋各水文要素的變化特性及各要素之間的相關聯(lián)系和影響機制。