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

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

1 研究背景

河流的水文要素變化過程是氣候變化和下墊面要素交互作用的結(jié)果，具有復雜的時空變化演變規(guī)律[1]。近年來，隨著全球極端水文、氣候事件頻發(fā)，了解河川徑流時空變化規(guī)律成為社會的迫切需求[2-3]。地球作為一個復雜的開放系統(tǒng)，太陽輻射是其重要的外部能量輸入和影響因素。研究表明，全球和區(qū)域降水在長時間尺度上受到太陽活動的影響，但其具體的影響過程和相互關(guān)系在不同區(qū)域有著不同的表現(xiàn)形式和規(guī)律特征[4]。Neff等[5]的研究論證了全新世以來，太陽活動與印度季風之間復雜的耦合關(guān)系。Zhang等[6]的研究也表明，太陽活動的變化對東亞區(qū)域降水和河流流量有著復雜的影響機制。竇睿音等[7]的研究表明，太陽黑子數(shù)與關(guān)中平原的整體干濕狀況的關(guān)系存在轉(zhuǎn)換現(xiàn)象，20世紀70年代之前二者為反向相關(guān)，之后二者之間的關(guān)系轉(zhuǎn)換為正相關(guān)。

河流流量的變化是區(qū)域氣象水文活動的綜合體現(xiàn)，反映了流域范圍內(nèi)整體的干濕狀況，因此流量變化成為氣候變化區(qū)域響應的指示器[8]。相關(guān)研究表明，盡管太陽黑子活動與河流流量之間存在相關(guān)性，但是由于氣象水文過程屬于非線性、非穩(wěn)定的復雜系統(tǒng)，如利用傳統(tǒng)統(tǒng)計方法分析其相關(guān)性可能會受到原始序列中噪音的干擾[9]，難以探究二者之間在不同時域和頻域上的復雜相關(guān)關(guān)系。因此，本文結(jié)合集合經(jīng)驗模態(tài)分解(Ensemble Empirical Mode Decomposition, EEMD)[10]和交叉小波變換(Cross Wavelet Transform, CWT)[11]，提出一種改進的相關(guān)分析方法——EEMD-CWT綜合分析，并將其應用于贛江長時間流量序列和同期太陽黑子數(shù)序列之間的相關(guān)性分析。文中結(jié)果一方面可以為贛江流域水資源和干濕狀況長期變化規(guī)律和影響因素的探討提供基礎資料，另一方面可以為提高氣象水文系統(tǒng)與外部強迫之間的非線性關(guān)系研究提供參考和思路。

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

2.1 流量數(shù)據(jù)與太陽黑子數(shù)據(jù)

贛江的流量數(shù)據(jù)來自1950-2015年南昌外洲水文站的實測資料。外洲水文站位于江西省南昌市范圍內(nèi)的贛江下游右岸(115.84°E、28.63°N)，控制面積為80 948 km2，占整個贛江流域范圍的99%以上[12]，因此外洲水文站的流量數(shù)據(jù)可以反映整個贛江流域的產(chǎn)流情況。外洲水文站設立于1949年10月，從1950年開始有連續(xù)的流量觀測數(shù)據(jù)，是目前贛江流域最長的流量數(shù)據(jù)序列(圖1(a))。

太陽黑子數(shù)據(jù)來自比利時皇家天文臺數(shù)據(jù)中心，該研究機構(gòu)是世界權(quán)威的太陽活動觀測中心，因此本文選取其發(fā)布的1950-2015年太陽黑子相對數(shù)(圖1(b))作為反映太陽活動狀況的序列(數(shù)據(jù)下載網(wǎng)址為http://www.sidc.be/silso/datafiles)。

圖1 1950-2015年贛江年均流量和太陽黑子數(shù)

由圖1可以觀察到，流量序列與太陽黑子序列之間的相關(guān)關(guān)系并不明顯。因此，利用SPSS 20軟件，對流量數(shù)據(jù)與太陽黑子數(shù)據(jù)進行初步的相關(guān)分析。結(jié)果表明，二者之間的Pearson相關(guān)系數(shù)為-0.118(P=0.347)。

2.2 研究方法

2.2.1EEMD方法 經(jīng)驗模態(tài)分解(EmpiricalModeDecomposition,EMD)是一種利用各信號周期不同的振幅頻率特征，將原始序列分解為一系列具有不同時間尺度本征模函數(shù)(IntrinsicModeFunction,IMF)的方法，被廣泛應用于水文、氣象、地球物理等領域的多尺度時頻分析[13-15]。集合經(jīng)驗模態(tài)分解(EEMD)則是EMD方法的改進算法，通過在原始信號中加入噪聲減弱了模態(tài)之間的混疊，從而提高了分解的有效性[16-17]。

2.2.2CWT方法 交叉小波變換(CrossWaveletTransform,CWT)是一種結(jié)合交叉譜分析和小波分析的新方法，可以同時在時域和頻域上分析兩個信號之間的相關(guān)性[11]。交叉小波能量譜可以反映出兩個信號經(jīng)過小波分析之后具有相同能量譜的區(qū)域，因此可以揭示出兩個信號在不同時頻域上相互作用的程度。其具體計算方法為：

(1)

2.2.3EEMD-CWT綜合分析 本文提出的EEMD-CWT綜合分析的思路是：首先利用EEMD方法將原始信號分解為不同尺度的互相正交的本征模函數(shù)(IMF)分量，以實現(xiàn)信號的平穩(wěn)化，降低原始信號中的噪音影響。然后利用多窗譜分析(Multi-TaperMethod,MTM)[18]確定原始信號中所存在的各種周期。最后，選擇具有相同物理意義的周期尺度下的本征模函數(shù)(IMF)分量進行交叉小波變換，得到二者的交叉小波能量譜和小波凝聚譜，并以此來分析兩個序列在不同時頻域上的相關(guān)性。

3 結(jié)果與分析

3.1 流量的多尺度特征和周期

對1950-2015年贛江流量序列進行EEMD分解(圖2)，共得到6個IMF分量，其中IMF1～5為不同時間尺度下流量序列的周期項，IMF6為整個序列的長期變化趨勢項。由圖2可以看出，IMF1分量的波動振幅最大，頻率也最高，但是周期性較差。這一信號代表了原始序列中的主要變化因素，是水文序列非線性和非平穩(wěn)特征的體現(xiàn)。EMD方法的提出者——Huang等[18]的研究指出，如果趨勢項為單調(diào)上升或單調(diào)下降函數(shù)則說明原始信號具有明顯的非平穩(wěn)性。由圖2可以清楚的觀察到IMF6為單調(diào)上升函數(shù)，這也論證了贛江近66年來的流量序列具有明顯的非平穩(wěn)性。IMF2～5這4個分量的周期性和平穩(wěn)性已經(jīng)得到了很大程度上的改善，表現(xiàn)出了比較明顯的周期變化規(guī)律。并且隨著分解尺度的遞進，其振動幅度和頻率逐漸下降，信號的平穩(wěn)性逐漸增強。如IMF4和IMF5已經(jīng)比較接近于正弦波。

(圖中縱坐標為相對值)

通過對1950-2015年贛江流量序列的多窗譜(MTM)分析表明(圖3)，近66年來贛江流量序列具有明顯的22、8、4.5和3.7a周期(通過α=0.05顯著性檢驗)，其中11a周期也比較明顯，但是未能通過相關(guān)檢驗。譜值最大的22a周期即為太陽活動的Hale周期[19]，也就是常見的太陽黑子11a周期的2倍。另外，準4～8a的周期很可能與ENSO的活動周期有關(guān)[20]。

(虛線表示95%置信區(qū)間)

3.2 太陽黑子的多尺度特征和周期

對1950-2015年太陽黑子數(shù)序列進行EEMD分解(圖4)，共得到6個IMF分量，其中IMF1～5為不同時間尺度下太陽黑子數(shù)序列的周期項，IMF6為整個序列的長期變化趨勢項。由圖4可以看出，太陽黑子數(shù)的變化過程較為平穩(wěn)，具有明顯的周期性特征。在IMF1分量中，除1950-1960年這段時間外，其余時間其振幅都較為平穩(wěn)。IMF2～5這4個分量具有非常明顯的周期變化特征，各個時間尺度下的變化過程都具有類似正弦波的結(jié)構(gòu)。IMF6趨勢項的變化過程顯示，1950-1970年為太陽黑子數(shù)上升階段，1970-2015年為下降階段。

(圖中縱坐標為相對值)

進一步的多窗譜分析(MTM)表明(圖5)，太陽黑子數(shù)序列具有顯著的36、11、5.4和3.2a周期(通過α=0.05顯著性檢驗)。其中太陽黑子數(shù)的11a周期表現(xiàn)最為明顯，其譜值明顯大于其他幾個周期。

(虛線表示95%置信區(qū)間)

3.3 流量與太陽黑子交叉小波分析

現(xiàn)代研究表明[21-22]，太陽活動異常對地球氣候系統(tǒng)、降水、流域流量的影響主要是通過影響環(huán)流和季風強度實現(xiàn)的。如過去120年以來南亞和東亞地區(qū)太陽輻射的增強通過影響沃克環(huán)流和哈德來環(huán)流，最終導致區(qū)域降水和流域流量的增加。因此可以通過太陽活動和流量之間的關(guān)系，進一步揭示二者在不同尺度下的相互影響過程。

通過前文的分析可知，贛江流量序列與太陽黑子數(shù)序列的相關(guān)系數(shù)僅為-0.118，且未能通過α=0.1的顯著性檢驗，不能直接利用未經(jīng)處理的原始序列進行二者之間的相關(guān)分析。但是，經(jīng)過EEMD分解后，各IMF分量的平穩(wěn)性有了很大程度的改善，周期性特征也更加明顯。與此同時，流量和太陽黑子的多窗譜分析(MTM)表明，二者都具有明顯的太陽活動周期，因此可以選取具有與太陽活動周期相關(guān)的IMF分量進行對比分析。以此為基礎，經(jīng)過交叉小波分析得到贛江流量和太陽黑子數(shù)之間的交叉小波能量譜和小波凝聚譜。

圖6為1950-2015年贛江流量序列IMF2分量與太陽黑子數(shù)序列的IMF2分量對比。由圖6可知，贛江流量變化與太陽黑子數(shù)序列之間的相關(guān)關(guān)系在不同時間段內(nèi)是不一樣的。其中，1950-1960年呈負相關(guān)、1961-1973年呈正相關(guān)、1974-1984年呈負相關(guān)、1985-1995年呈正相關(guān)、1996-2007年為負相關(guān)(部分時段相差1/4個周期)、2008-2015年為正相關(guān)(部分時段相差1/4個周期)。近66年來，大體上二者以11a為周期，在負相關(guān)和正相關(guān)之間交替變換。

(圖中縱坐標為相對值)

圖7為1950-2015年贛江流量與太陽黑子數(shù)序列的交叉小波譜和小波凝聚譜，交叉小波能量譜中顏色越偏紅表示能量譜密度值越大。由圖7(a)可以觀察到在11a周期附近贛江流量與太陽黑子數(shù)的相關(guān)關(guān)系最為顯著，其中黑色實線所圍合區(qū)域表示通過了95%置信度檢驗，灰色陰影區(qū)則表示可能會受到邊界效應的影響。圖中→表示二者之間變化的位相一致(正相關(guān))，←表示位相相反(負相關(guān))，↑和↓則表示二者之間提前或落后1/4個周期。在圖7(a)中的11a周期尺度上，1950-2015年二者之間由負相關(guān)逐漸順時針轉(zhuǎn)為相差1/4個周期，再轉(zhuǎn)為正相關(guān)。小波凝聚譜(圖7(b))中相關(guān)系數(shù)較高的區(qū)域與交叉小波能量譜(圖7(a))中基本一致，也是在11a尺度上達到最大值。其中1950-2000年，二者的相關(guān)系數(shù)在0.6以上，2001-2015年二者的相關(guān)系數(shù)也在0.4以上。小波凝聚譜可以發(fā)現(xiàn)在32a周期尺度上，二者的相關(guān)系數(shù)值在整個研究時段都在0.8以上，但是考慮到本研究的長度只有66a，而且由于邊界效應的存在，這一相關(guān)關(guān)系可能并不具有實際意義。

圖7 1950-2015年贛江流量與太陽黑子數(shù)序列的交叉小波譜和小波凝聚譜

4 結(jié) 論

本文基于贛江外洲水文站1950-2015年的長時間流量序列和同期太陽黑子數(shù)序列，結(jié)合EEMD及CWT綜合分析方法對流量與太陽活動之間的多尺度相關(guān)關(guān)系進行了分析。研究結(jié)果表明：

(1)EEMD-CWT綜合分析可以有效利用EEMD方法將原始信號進行分解，得到具有平穩(wěn)性的IMF分量，并借助多窗譜分析(MTM)識別出信號中的不同周期，選擇具有相同物理意義的贛江流量序列IMF分量和太陽黑子IMF分量，再利用交叉小波變換(CWT)得到二者在時域和頻域上的相關(guān)關(guān)系變化過程。

(2)1950-2015年贛江流量序列的EEMD分解結(jié)果表明，該序列具有明顯的非平穩(wěn)性，其中IMF1分量的波動振幅最大，頻率最高，代表了原始序列中的主要變化因素，是水文序列非線性和非平穩(wěn)特征的體現(xiàn)。

(3)1950-2015年太陽黑子數(shù)序列的EEMD分解結(jié)果表明，太陽黑子數(shù)的變化過程較為平穩(wěn)，具有明顯的周期性特征。趨勢項的變化過程顯示，1950-1970年為太陽黑子數(shù)上升階段，1970-2015年為下降階段。

(4)多窗譜分析(MTM)表明，近66年來贛江流量序列具有22、8、4.5和3.7a周期，太陽黑子數(shù)序列具有36、11、5.4和3.2a周期。

(5)交叉小波分析(CWT)表明，在11a周期尺度上，1950-2015年贛江流量與太陽黑子數(shù)之間由負相關(guān)逐漸順時針轉(zhuǎn)為相差1/4個周期，再轉(zhuǎn)為正相關(guān)。其中1950-2000年，二者的相關(guān)系數(shù)在0.6以上，2001-2015年相關(guān)系數(shù)在0.4以上。

[1]GATESWL.Anoverviewoftheresultsoftheatmosphericmodelintercomparisonproject(AMIPI)[J].BulletinoftheAmericanMeteorologicalSociety,1999，80(2)：29-55.

[2] 夏 軍，談 戈. 全球變化與水文科學新的進展與挑戰(zhàn)[J]. 資源科學,2002，24(3)：1-7.

[3] 范垂仁，夏 軍，張利平，等. 中國水旱災害長期預報[M]. 北京：中國水利水電出版社，2008.

[4] 張?zhí)m生，方修琦，任國玉. 全球變化[M]. 北京：高等教育出版社，2000.

[5]NEFFU，BURNSSJ，MANGINIA，etal.Strongcoherencebetweensolarvariabilityandthemonsooninomanbetween9and6KyrAgo[J].Nature,2001，411(6835)：290-293.

[6]ZHANGP，CHENGH，EDWARDSRL，etal.Atestofclimate，sun，andculturerelationshipsfroman1810-Yearchinesecaverecord[J].Science,2008，322(5903)：940-942.

[7] 竇睿音，延軍平. 關(guān)中平原太陽黑子活動周期與旱澇災害的相關(guān)性分析[J]. 干旱區(qū)資源與環(huán)境,2013，27(8)：76-82.

[8]CHENF，YUANYJ，WEIWS，etal.Tree-ringrecordedhydroclimaticchangeinTienshanMountainsduringthepast500years[J].QuaternaryInternational,2015,358：35-41.

[9]YEHJ，SHIEHJ，HUANGNE.Complementaryensembleempiricalmodedecomposition：anovelnoiseenhanceddataanalysismethod[J].AdvancesinAdaptiveDataAnalysis,2010，2(2)：135-156.

[10]WUZhaohua，HUANGNE.EnsembleEmpiricalModeDecomposition：Anoise-assisteddataanalysismethod[J].AdvancesinAdaptiveDataAnalysis,2009,1(1)：1-41.

[11]GRINSTEDA，MOOREJC，JEVREJEVAS.Applicationofthecrosswavelettransformandwaveletcoherencetogeophysicaltimeseries[J].NonlinearProcessesinGeophysics,2004，11(5/6)：561-566.

[12] 顧朝軍，穆興民，高 鵬，等. 贛江流域徑流量和輸沙量的變化過程及其對人類活動的響應[J]. 泥沙研究,2016,41(3)：38-44.

[13]JIF，WUZ，HUANGJ，etal.Evolutionoflandsurfaceairtemperaturetrend[J].NatureClimateChange,2014，4(6)：462-466.

[14]CHENYangkang，MAJitao.RandomnoiseattenuationbyF-Xempiricalmodedecompositionpredictivefiltering[J].Geophysics,2013，79(3)：81-91.

[15]BREAKERLC，LOORHR，CARROLLD.Trendsinseasurfacetemperatureoffthecoastofecuadorandthemajorprocessesthatcontributetothem[J].JournalofMarineSystems,2016，164：151-164.

[16]WUZ，HUANGNE，CHENX.Themulti-dimensionalensembleempiricalmodedecompositionmethod[J].AdvancesinAdaptiveDataAnalysi,2009，1(3)：339-372.

[17]MANNME，LEESJM.Robustestimationofbackgroundnoiseandsignaldetectioninclimatictimeseries[J].ClimaticChange,1996，33(3)：409-445.

[18]HUANGNE，SHENZ，LONGSR，etal.Theempiricalmodedecompositionandthehilbertspectrumfornonlinearandnon-stationarytimeseriesanalysis[J].ProceedingsoftheRoyalSocietyofLondonA(Mathematical，PhysicalandEngineeringSciences),1998，454(1971)：903-995.

[19]RINDD.Thesun'sroleinclimatevariations[J].Science,2002，296(5568)：673-677.

[20] 王鐘睿，錢永甫. 江淮梅雨的多尺度特征及其與厄爾尼諾和大氣環(huán)流的聯(lián)系[J]. 大氣科學學報. 2004，27(3)：317-325.

[21] 葛全勝，劉路路，鄭景云，等. 過去千年太陽活動異常期的中國東部旱澇格局[J]. 地理學報，2016，71(5):707-717.

[22]BHATTACHARYYAS，NARASIMHAR.RegionaldifferentiationinmultidecadalconnectionsbetweenIndianmonsoonrainfallandsolaractivity[J].JournalofGeophysicalResearchAtmospheres，2007，112(24):177-180.