萬(wàn)智巍，賈玉連，洪祎君，蔣梅鑫

(1.江西師范大學(xué)地理與環(huán)境學(xué)院鄱陽(yáng)湖濕地與流域研究教育部重點(diǎn)實(shí)驗(yàn)室，江西南昌330022；2.中國(guó)科學(xué)院地理科學(xué)與資源研究所陸地表層格局與模擬院重點(diǎn)實(shí)驗(yàn)室，北京100101)

0 引 言

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

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

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

1 數(shù)據(jù)與方法

1.1 徑流數(shù)據(jù)

贛江的徑流數(shù)據(jù)來(lái)自1950年～2015年南昌外洲水文站的實(shí)測(cè)數(shù)據(jù)。該水文站位于江西省南昌市境內(nèi)的贛江下游右岸，集水面積為80 948 km2，占整個(gè)贛江流域的99.6%[17]。外洲站設(shè)立于1949年10月，自1950年開(kāi)始有連續(xù)的水文觀測(cè)記錄，是目前贛江流域具有最長(zhǎng)連續(xù)水文觀測(cè)資料的站點(diǎn)。

1.2 研究方法

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

(1)

式中，cij(t)為第i次加入白噪聲后得到的第j個(gè)分量，也稱(chēng)之為IMF(本征模函數(shù))；rij(t)為第j次分解信號(hào)后得到的剩余趨勢(shì)項(xiàng)。

利用線性?xún)A向和MK趨勢(shì)檢驗(yàn)對(duì)徑流序列進(jìn)行判斷。線性?xún)A向估計(jì)的原理是利用一元線性回歸方程擬合時(shí)間序列和觀測(cè)值，并以回歸系數(shù)作為判斷序列長(zhǎng)期趨勢(shì)的方法[18]。其具體求解公式為

y=ax+b

(2)

(3)

式中，y為徑流；x為年份；a為回歸系數(shù)(線性?xún)A向估計(jì)值)。

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

2 贛江徑流變化特征

2.1 年際變化與趨勢(shì)

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

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

由圖1可以發(fā)現(xiàn)，贛江年徑流變化波動(dòng)范圍較大，最小值出現(xiàn)于1963年，徑流量為236.07億m3；最大值出現(xiàn)于1973年，為1147.87億m3；其多年平均值為682.46億m3。線性?xún)A向估計(jì)的結(jié)果表明，贛江徑流變化具有上升趨勢(shì)，其上升速率為6.5×108m3/10 a。MK趨勢(shì)檢驗(yàn)的結(jié)果表明，其Z統(tǒng)計(jì)量為0.653。由于其Z值<1.64，因此未能通過(guò)相關(guān)顯著性檢驗(yàn)，也就是說(shuō)從長(zhǎng)期來(lái)看贛江徑流變化的上升趨勢(shì)并不顯著。

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

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

圖3 1950年～2015年徑流序列MK突變檢驗(yàn)

2.2 多尺度變化特征

2.2.1 模態(tài)特征

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

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

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

表1 IMF分量和趨勢(shì)分量的主要統(tǒng)計(jì)值

2.2.2 周期特征

為了進(jìn)一步揭示出不同時(shí)間尺度下贛江徑流變化的周期特征，同時(shí)結(jié)合上一節(jié)分析結(jié)果中所得出的CEEMD分解前4項(xiàng)方差之和占到序列總方差的99.5%這一特征，利用功率譜分析方法對(duì)IMF1～I(xiàn)MF4進(jìn)行周期分析，以確定其所包含的全部周期信息。結(jié)果顯示(見(jiàn)圖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%置信區(qū)間)

3 討 論

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

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

4 結(jié) 論

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

(1)線性?xún)A向估計(jì)的結(jié)果表明，贛江徑流變化上升速率為6.5×108m3/10 a。MK趨勢(shì)檢驗(yàn)的結(jié)果表明，其Z統(tǒng)計(jì)量為0.653，上升趨勢(shì)并不顯著。

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

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

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