段志鵬，李繼清

(華北電力大學可再生能源學院，北京 102206)

0　引　言

北京市是我國水資源緊缺最嚴重的城市之一。在全球氣候變暖的背景下，北京市在21世紀以來，全年降水明顯減少，更是加劇了北京市的缺水危機，因此研究北京市降水的變化規(guī)律以及未來的演變趨勢，成為人們普遍關(guān)注的中心問題。目前已經(jīng)有許多學者對北京市降水的變化規(guī)律及演變趨勢等方面做了許多研究，朱龍騰[1]等采用經(jīng)驗模態(tài)分解(EMD)方法分析了北京市1951～2009年降水時間序列的多尺度振蕩變化，分析結(jié)果表示短期內(nèi)北京市的年降水量將在波動中不斷減少；孫振華[2]等通過對北京市16個雨量站55年的降水資料進行統(tǒng)計計算和頻率分析，探討了北京市降水的時空分布特征及變化趨勢，發(fā)現(xiàn)20世紀90年代之后北京市連續(xù)干旱，2008年之后年降水可能會進入一個相對較豐的時期；王文[3]等利用1961-2008年北京11個臺站的氣候觀測資料，分析了過去48 a北京城區(qū)和郊區(qū)的氣溫和降水變化趨勢，發(fā)現(xiàn)無論郊區(qū)和城區(qū)的年降水量都是呈減少趨勢。上述研究成果從不同角度分析了北京市氣溫、降水的周期和變化特征，不僅有利于深入理解氣候變化過程與規(guī)律，進行氣候預測，更對北京市社會經(jīng)濟發(fā)展具有重要意義。

氣候系統(tǒng)是非線性復雜系統(tǒng)，相關(guān)研究表明，氣候系統(tǒng)除了具有非線性和非平穩(wěn)性，還有層次性，許多大小不一的時空尺度構(gòu)成了多層次結(jié)構(gòu)[4]，傳統(tǒng)趨勢分析中，大尺度循環(huán)和趨勢變化很可能混合在一起，不能分辨出趨勢變化還是周期震蕩。極點對稱模態(tài)分解(ESMD)方法是近年來發(fā)展的局部自適應時間序列分析技術(shù)，是在經(jīng)驗模態(tài)分解方法基礎上改進的一種新的非線性、非平穩(wěn)時間序列分析方法。較之常用的小波分析方法，ESMD可以根據(jù)時間序列局部時變特征進行自適應地時頻分解，完全擺脫了傅里葉變換的束縛，并可得到極高的時頻分辨率，非常適合對非平穩(wěn)、非線性時間序列進行分析。與廣泛應用于氣候水文要素分析的EMD方法相比，該方法解決了EMD方法中的“模態(tài)混疊”(或稱“頻率交叉”)問題，能夠提供更好的尺度分解，模態(tài)分量可以準確體現(xiàn)時間序列的波動特征和趨勢變化，在氣候數(shù)據(jù)分析方面具有較大優(yōu)勢。

鑒于此本文采用ESMD分解對北京市1951-2015年的年降水序列進行多時間尺度和突變、趨勢分析，并結(jié)合小波分析、Mann-Kendall 法和重標極差(R/S)分析法對ESMD方法的結(jié)果進行對比驗證，為今后北京市降水預報和水資源規(guī)劃管理提供參考。

1　研究方法

1.1　經(jīng)驗模態(tài)分解的基本原理

1998年，黃鍔提出了經(jīng)驗模態(tài)分解(Empirical Mode Decomposition，EMD)方法，這種新的數(shù)據(jù)處理方法能夠在時域和頻域上提取數(shù)據(jù)的信息，很好的保留數(shù)據(jù)本身的特性，具有良好的自適應性，非常適合處理非線性、非平穩(wěn)信號[5]。EMD方法是把不同周期的波動從原信號中分離出來，并且該波動是平穩(wěn)的，最后得到的是趨勢分量。不同尺度的波動被定義成為本征模態(tài)函數(shù)(Intrinsic Mode Function，IMF)。經(jīng)驗模態(tài)分解的計算過程如下：

(2)計算h1(t)=x(t)-m1(t)；

(3)將h1(t)作為新的信號序列，重復上述操作直到h1(t)成為一個零均值過程；

(4)得出零均值過程h1(t)后，將其作為第一個IMF分量c1(t)，它表示原始信號最高頻的分量；

(5)從原始信號x(t)中減去c1(t)，得到去除一個高頻分量的新信號過程r1(t)，即：

r1(t)=x(t)-c1(t)

(6)將r1(t)作為原始數(shù)據(jù)，重復(1)～(4)步驟，依次可得到c2(t)、c3(t)、c4(t)等，一直到滿足迭代條件為止，則原始信號重構(gòu)為：

(1)

式中：rn(t)為趨勢項，代表信號的平均趨勢。

1.2　極點對稱模態(tài)分解原理

極點對稱模態(tài)分解(Extreme-point Symmetric ModeDecomposition，ESMD)方法是經(jīng)驗模態(tài)分解方法的新發(fā)展，該方法解決了EMD方法中篩選終止和“模態(tài)混疊”(或稱“頻率交叉”)問題，并在氣候數(shù)據(jù)分析方面具有較大優(yōu)勢[6]：①擅長尋找變化趨勢，能夠觀測序列中分離出年際變化趨勢和氣候變化總趨勢，有助于探究全球氣候變暖問題；②擅長異常診斷，能夠從分解模態(tài)中發(fā)現(xiàn)異常時段與頻段，有利于氣候異常研究；③擅長時-頻分析，先進的直接插值法能夠直觀地分析各時間尺度上的頻率變化。

ESMD方法是由兩部分組成：第一部分是模態(tài)分解，可產(chǎn)生數(shù)個模態(tài)與一條最佳自適應全局均線；第二部分是時-頻分析，利用直接插值法計算固有模態(tài)的瞬時頻率，分析各時間尺度上的頻率變化，判斷各時間尺度下突變發(fā)生的時間。ESMD方法模態(tài)分解過程如下：

(1)找出序列X的所有極大值和極小值點，記作Ei(1≤i≤n)。

(2)用線段連接所有相鄰的極點并將線段中點依次記為Fi(1≤i≤n-1)，并通過一定方式在左右兩端添補邊界中點F0和Fn。

(3)利用n+1個中點構(gòu)建p條內(nèi)插曲線L1，…，Lp(p≥1)，計算其均值曲線L*=(L1+…+Lp)/p。

(4)對X-L*序列重復上述步驟，直到|L*|≤ε(ε是允許誤差)或者篩選次數(shù)達到預設的最大值K，得到第一個經(jīng)驗模M1。

(5)對剩余序列X-M1重復上述4個步驟，直到剩余序列R只剩于一定數(shù)量的極點，便可分別得到經(jīng)驗模M2，M3,…。

(6)在限定區(qū)間[Kmin,Kmax]內(nèi)改變最大篩選次數(shù)K值，重復上述5個步驟。然后計算方差比率σ/σ0，并畫出它隨K的變化圖，在圖中找出σ/σ0最小值對應的K0，以K0作為限制條件再次重復上述5個步驟，最后剩余項R就是序列X的自適應全局均線。

經(jīng)過分解，原始的時間序列X可表示為X=∑Mi+R，即利用ESMD將時間序列X分解成一系列經(jīng)驗模態(tài)和一個剩余變量。

在第二部分時-頻分析中，考慮到原有的希爾伯特變換把原本離散的信號轉(zhuǎn)換成解析函數(shù)處理，這樣的處理方式不可避免的要受到各種數(shù)學理論的限制。面對離散數(shù)據(jù)，在分析過程中應當尊重其離散性特征，故此提出了針對數(shù)據(jù)的“直接插值法”，借此不但可以直觀地體現(xiàn)各模態(tài)的振幅與頻率的時變性，還可明確地獲知總能量變化[7]。直接插值法基本思路如下：

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

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

(3)將這些局部周期值取倒數(shù)得到局部頻率，再做三次樣條插值得到光滑的時間-頻率變化曲線。

2　應用實例

2.1　研究區(qū)概況

北京市地處華北平原西北端，海河流域中部，土地面積16 410 km2，山區(qū)面積約占總面積的61%。地理環(huán)境特殊，海拔高度變差大，山地和平原之間過渡急劇，界限清晰。北京市屬于典型的北溫帶半濕潤大陸性季風氣候，夏季高溫多雨，冬季寒冷干燥，降水季節(jié)分配不均。地形的復雜多變造成了北京地區(qū)天氣狀況的多樣性。暴雨和強對流天氣是北京地區(qū)夏季(6-8月)主要的災害性天氣[8,9]，導致北京地區(qū)旱澇災害頻繁。如圖1所示為北京市1951-2015近65 a的降水過程。其中，年最大降水量為1 407 mm，發(fā)生在1957年；年最小降水量為267 mm，出現(xiàn)在1965年；多年平均降水量為589.75 mm。對北京市的年降水量進行線性擬合，可以看出從1951-2015年年平均降水量整體呈下降趨勢，線性降水減少速率為3.47 mm/a。

圖1　北京市1951-2015年降水時間序列Fig.1 Annual precipitation time series from 1951 to 2015 at Beijing station

本文研究所用年降雨資料來源于中國氣象局國家氣候中心提供的北京市1951-2015年月降雨資料，是經(jīng)過該機構(gòu)整編審查后發(fā)布的，故資料的可靠性基本可以得到保證。對資料的代表性進行分析，如圖2所示，繪制北京市年降雨量模比系數(shù)差積曲線，可以看出隨時間的增長，曲線變幅越來越小，年降雨系列在50年以上時，模比系數(shù)差積曲線穩(wěn)定趨近于“1”，說明此系列資料具有一定的穩(wěn)定性。

圖2　北京市年降雨量模比系數(shù)差積曲線Fig.2 Beijing annual precipitation model ratio coefficient curve

2.2　降水變化的多尺度分析

ESMD方法具有自適應性和基于信號的局部變化特性，適用于非線性、非平穩(wěn)序列的時頻分析，與其他方法相比可以更好地提取非線性序列的變化趨勢。因此利用ESMD方法對北京市1951-2015年年平均降水量的時間序列進行了分解，計算出當方差比率最小時最佳篩選次數(shù)為26次，此時得到的分解結(jié)果最好，為4個模態(tài)分量(模態(tài)1～4)以及一個趨勢余量R，如圖3所示。

圖3　降水時間序列ESMD分解的模態(tài)分量及趨勢項Fig.3 Modal components and trend of precipitation time series by ESMD

為驗證分解結(jié)果的可靠性，用分解得到的模態(tài)1～4和趨勢余量R進行合成，可得到一個重構(gòu)序列，發(fā)現(xiàn)重構(gòu)序列與原始的降水序列完全吻合，這說明ESMD方法分解結(jié)果是可信的。圖2中ESMD分解的4個模態(tài)分量各自都具有獨立的代表性，不會在同一時刻出現(xiàn)時間尺度相同的振蕩，且每個模態(tài)分量一般都具有物理意義，它們各自反映了原序列中固有的不同特征尺度的振蕩。為解析降水序列固有的不同特征尺度的振蕩，利用周期圖法[10]估算各個分解模態(tài)的平均周期，進而將每種尺度信號波動頻率和振幅對原數(shù)據(jù)總體特征影響用方差貢獻率表示出來，如表1所示。

表1　降水序列各分量的方差貢獻率Tab.1 Contribution rates of ESMD decomposition for mean precipitation

表1給出了不同尺度模態(tài)分量所表征的不同時間尺度波動的平均周期和各模態(tài)分量的方差貢獻率。由降水序列的各模態(tài)分量可以看出，在年際尺度上，北京市年平均降水以2.6 a(模態(tài)1)和4.3 a(模態(tài)2)的周期振蕩為主；在年代際尺度上，則以14 a(模態(tài)3)和21.7 a(模態(tài)4)的周期振蕩為主。

結(jié)合表1和圖2可以看出：北京市65 a的降水序列包含多個時間尺度特征。其中模態(tài)1[圖2(a)]表示的2.6 a周期方差貢獻率最大，達到了32.47%，信號振蕩非常明顯，表明降水量變化呈現(xiàn)減少—增加—減少的循環(huán)變化，基本反映了20世紀70年代以前北京市降水豐枯交替的變化情況；模態(tài)2[圖2(b)]表示的4.3 a周期方差貢獻率約為16.74%，20世紀70年代以前振幅波動較大，而70年代以后振幅逐漸減小，趨于平穩(wěn)，與模態(tài)1表現(xiàn)出相似特征；模態(tài)3[圖2(c)]表示的14 a周期方差貢獻率約為27.5%，僅次于模態(tài)1，振幅在20世紀60年代末期至20世紀90年代中期變化平穩(wěn)，波動相對較小，而在20世紀50年代初期至60年代末期和21世紀之后，振幅明顯增大，說明相比于其他時期，這兩個時期的降水量變化幅度較大；模態(tài)4[圖2(d)]表示的21.7年周期方差貢獻率僅為1.37%，整個時間跨度上波動幅度比較穩(wěn)定。對比模態(tài)3和模態(tài)4可以發(fā)現(xiàn)，模態(tài)3的振幅變化大而模態(tài)4的振幅變化顯著變小，這表明降水異常主要發(fā)生在周期為準14 a的時間尺度上。趨勢項R[圖2(e)]為最佳自適應全局均線對應于年際降水變化，其方差貢獻率為21.92%，反映了過去65 a降水序列的整體變化趨勢，從圖2(e)中可以看出，趨勢項在整個時間尺度上呈現(xiàn)非線性的逐漸下降趨勢。1951-1970年降水量逐年減少，1970-1990年降水量總體變化不大，1990-2015年降水量又逐年下降，且下降趨勢較為明顯。

小波分析是近年來十分熱門的一個前沿領(lǐng)域，因其對信號處理具有特殊優(yōu)勢，已被廣泛應用于氣象和氣候序列的時頻結(jié)構(gòu)分析中[11]。為了對比ESMD方法和小波分析在信號處理中的性能，采用水文氣象上常用的Morlet小波變換對北京市的年降水距平序列進行分析，并且為了減弱序列開始點和結(jié)束點附近的邊界效應，采用對稱外延法將降水序列擴展[12]。如圖4所示為北京市近65 a降水距平序列的小波系數(shù)實部時頻分布圖。該圖反應了降水序列各種周期的強弱和分布，小波系數(shù)為正時，表示降水量相對偏多，值越大降水越多；小波系數(shù)為負時表示降水量偏少。從圖3中可以看出北京市降水存在多時間尺度特征，降水量變化存在準2、5～8、10～15和準25 a的周期。整個時間跨度上正負相位交替出現(xiàn)，說明降水量經(jīng)歷豐枯交替變化。綜上所述，表2給出了小波與ESMD分析周期結(jié)果對比，二者結(jié)果有一定差異，而且呈現(xiàn)出周期越長差異越大的規(guī)律。究其原因，小波分析是以傅里葉變換作為理論依據(jù)，存在一定的局限性，如小波基的選取、固定的基函數(shù)和恒定的多分辨率等問題[13]。與小波分析相比，ESMD方法擺脫了傅里葉變換理論的束縛，具有較強的靈活性和自適應性，依據(jù)數(shù)據(jù)自身特點進行分解，更有助于探索事物的內(nèi)在規(guī)律。

圖4　小波系數(shù)實部時頻分布圖Fig.4 The real part distribution of wavelet coefficients

表2　ESMD與小波分析周期結(jié)果對比　aTab.2 Comparison of cycle results by ESMD and Wavelet Analysis

2.3　不同時間尺度下降水突變分析

將分解得到的4個固有模態(tài)利用直接插值法繪出其瞬時頻率與瞬時振幅的時變圖，如圖5所示，圖5中F代表頻率，A代表振幅。這樣的分布圖可以明確表達出模態(tài)的振蕩強度與變化快慢程度。觀察F1和A1可知，在2.6 a尺度下，突變點較多，降雨偏多和偏少交替相對比較頻繁，如1957年時發(fā)生了突然的低頻、大振幅振蕩，對比原始數(shù)據(jù)可以發(fā)現(xiàn)，1957年北京市降水量僅487 mm，而1956年降水量達到了1 116 mm，下降了約56.4%；對比F2和A2，4.3 a尺度下，突變點在1961年和1995年；而且F1和F2存在一個頻率恒定的時期，這是其他方法不可見的結(jié)果。F3和A3對比可知，在14 a尺度下，在1988年附近發(fā)生突變；觀察F4和A4，發(fā)現(xiàn)在21.7 a的年代際尺度上，20世紀70年代以后的模態(tài)頻率持續(xù)增高，而振幅卻一直減小，說明在21.7年時間尺度上北京地區(qū)年降雨量持續(xù)減少。

圖5　降水數(shù)據(jù)各模態(tài)頻率(F)與振幅(A)的時變圖Fig.5 Time-varying figure for each modal frequency and amplitude of precipitation data

Mann-Kendall 法(簡稱M-K) 是一種非參數(shù)統(tǒng)計檢驗方法，M-K法能很好地揭示時間序列的突變及其趨勢變化特征[14]。為了驗證北京市年降水量突變的顯著性，作M-K檢驗曲線，可以確定突變發(fā)生在1990年附近，突變發(fā)生后北京市降水量呈現(xiàn)減少的趨勢。對比兩種方法的突變分析結(jié)果，不難看出ESMD方法具有更大的優(yōu)勢。M-K方法僅能找到整個時間序列的突變年份，而ESMD方法時頻分析通過不同模態(tài)瞬時頻率與瞬時振幅的對應關(guān)系圖，可以找出不同時間尺度下的突變年份，分析結(jié)果更加詳細全面，特別對局部時間序列的降水變化有較好分析，更好地為水文預報提供參考。

2.4　降水變化的趨勢分析

圖6為運用ESMD方法對降水序列進行不同時間尺度的分解后趨勢項對原始數(shù)據(jù)擬合情況，表示降水序列整體的變化趨勢。從圖6中可以看出，在整個時間尺度上呈現(xiàn)非線性的逐漸下降趨勢，預測到北京市降水量短期內(nèi)呈現(xiàn)波動中下降的趨勢。

圖6　最佳自適應全局均線對降雨數(shù)據(jù)的擬合情況圖Fig.6 The fitting of the best adaptive global mean curve to the rainfall data

R/S分析法最初由英國水文學家赫斯特(Hurst,1951年)在研究尼羅河水壩工程時提出，該方法認為降水量的變化是具有自相似性和長程依賴性的時間序列，Hurst指數(shù)H是定量描述時間序列自相似性與長程依賴性的有效方法[15,16]。利用R/S分析法側(cè)面驗證ESMD方法對于趨勢判斷的正確性。北京市年降水的Hurst指數(shù)H=0.649 38，大于0.5，說明年降水序列具有持續(xù)性，預測未來降水將持續(xù)減少。因此運用ESMD方法對降水序列的趨勢判斷具有可靠性，這種趨勢變化特征可以為北京市水資源的合理利用和有效管理提供科學的依據(jù)。

3　結(jié)　論

ESMD是一種適用于非線性、非平穩(wěn)序列的信號分析方法，將其應用于氣候要素時間序列，可提取氣候變化的固有時間尺度，不但能夠從數(shù)年的觀測序列中分離出年際和年代際變化趨勢，還能夠從其中分離出氣候變化的總趨勢。本文利用ESMD方法對北京市1951-2015年降水序列進行分解，可以得到以下結(jié)論：

(1)北京市65 a的降水序列包含多個時間尺度特征。在年際尺度上具有2.6 a和4.3 a的周期變化，在年代際尺度上具有14 a和21.7 a的周期變化，其中年際周期變化在北京市整體降水量變化中占據(jù)主導地位。

(2)不同時間尺度下突變發(fā)生的時間不同。在2.6 a尺度下，突變點較多，降雨偏多和偏少交替相對比較頻繁；4.3 a尺度下，突變發(fā)生在1961年和1995年附近；14 a尺度下，在1988年附近發(fā)生突變。

(3)北京市降水量在研究時期內(nèi)呈現(xiàn)“下降—平穩(wěn)—下降”的非線性下降趨勢，且后期下降較為明顯，因此預測未來幾年內(nèi)北京市年降水量仍將持續(xù)減少。

在全球變暖背景下，自然因素和人類活動影響相互疊加，水文時間序列的變化規(guī)律更加復雜，本文通過對北京市降水序列的研究表明ESMD方法在水文序列的多時間尺度和突變、趨勢分析上具有較好的效果。但是目前的研究僅處于初步階段，如何確定各個模態(tài)分量的影響因子，計算其權(quán)重并最終應用到水文時間序列預測中，將是以后研究的重點工作。

□

參考文獻：

[1]朱龍騰，陳遠生，李璐，等．1951-2009年北京市降水變化情勢分析[J]．水資源保護，2012，28(3)：42-46．

[2]孫振華，馮紹元，楊忠山，等．1950-2005年北京市降水特征初步分析[J]．灌溉排水學報，2007，(2)：12-16．

[3]王文，張薇，蔡曉軍.近50 a來北京市氣溫和降水的變化[J]．干旱氣象，2009，27(4)：350-353．

[4]薛春芳，侯威，趙俊虎，等．集合經(jīng)驗模態(tài)分解在區(qū)域降水變化多尺度分析及氣候變化響應研究中的應用[J]．物理學報，2013，62(10)：109-203．

[5]Huang N E，Shen Z，Long S R，et al.．The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis [J]．Proceedings：Mathematical Physical and Engineering Sciences，1998，454:903-995．

[6]王金良，李宗軍．可用于氣候數(shù)據(jù)分析的ESMD方法[J]．氣候變化研究快報，2014，(3):1-5．

[7]王金良，李宗軍．極點對稱模態(tài)分解方法——數(shù)據(jù)分析與科學探索的新途徑[M]．北京：高等教育出版社，2015：7-16．

[8]北京志編委．北京志地質(zhì)礦產(chǎn)水利氣象卷水利志[M]．北京：北京出版社，2000．

[9]北京市水利局．北京水旱災害[M]．北京：中國水利水電出版社，1999．

[10]李強，吳健，許正文，等．利用EMD方法提取太陽活動周期成分[J]．空間科學學報，2007，27(1):1-6 ．

[11]魏鳳英．現(xiàn)代氣候統(tǒng)計診斷與預測技術(shù)[M]．2版.北京：氣象出版社，2007:36-104．

[12]王文圣，丁晶，向紅蓮．小波分析在水文學中的應用研究及展望[J]．水科學進展，2002，13(4)：515-520．

[13]楊永峰，吳亞峰．經(jīng)驗模態(tài)分解在振動分析中的應用[M]．北京：國防工業(yè)出版社，2013:52-53．

[14]白晶，延軍平，蘇坤慧．1958- 2007 年秦嶺南北氣候變化的差異性分析[J]．陜西師范大學學報， 2010，38(6): 98-105．

[15]常遠勇，侯西勇，毋亭，等．1998-2010年全球低緯度降水時空特征分析[J]．水科學進展，2012，(4)．

[16]江田漢，鄧蓮堂．Hurst指數(shù)估計中存在的若干問題：以在氣候變化研究中的應用為例[J]．地理科學，2004，24(2)：177-182．