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

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

1 研究背景

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

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

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

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

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

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

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

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

2.2 研究方法

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

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

(1)

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

3 結(jié)果與分析

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

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

(圖中縱坐標(biāo)為相對(duì)值)

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

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

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

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

(圖中縱坐標(biāo)為相對(duì)值)

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

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

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

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

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

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

(圖中縱坐標(biāo)為相對(duì)值)

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

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

4 結(jié) 論

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

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

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

(3)1950-2015年太陽黑子數(shù)序列的EEMD分解結(jié)果表明，太陽黑子數(shù)的變化過程較為平穩(wěn)，具有明顯的周期性特征。趨勢項(xiàng)的變化過程顯示，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ù)之間由負(fù)相關(guān)逐漸順時(shí)針轉(zhuǎn)為相差1/4個(gè)周期，再轉(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] 夏 軍，談 戈. 全球變化與水文科學(xué)新的進(jìn)展與挑戰(zhàn)[J]. 資源科學(xué),2002，24(3)：1-7.

[3] 范垂仁，夏 軍，張利平，等. 中國水旱災(zāi)害長期預(yù)報(bào)[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)中平原太陽黑子活動(dòng)周期與旱澇災(zāi)害的相關(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] 顧朝軍，穆興民，高 鵬，等. 贛江流域徑流量和輸沙量的變化過程及其對(duì)人類活動(dòng)的響應(yīng)[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]. 大氣科學(xué)學(xué)報(bào). 2004，27(3)：317-325.

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

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