基于極點對稱模態分解的北京市降水特征分析

段志鵬，李繼清

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

0　引　言

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

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

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

1　研究方法

1.1　經驗模態分解的基本原理

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

(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)作為原始數據，重復(1)～(4)步驟，依次可得到c2(t)、c3(t)、c4(t)等，一直到滿足迭代條件為止，則原始信號重構為：

(1)

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

1.2　極點對稱模態分解原理

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

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

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

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

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

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

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

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

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

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

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

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

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

2　應用實例

2.1　研究區概況

北京市地處華北平原西北端，海河流域中部，土地面積16 410 km2，山區面積約占總面積的61%。地理環境特殊，海拔高度變差大，山地和平原之間過渡急劇，界限清晰。北京市屬于典型的北溫帶半濕潤大陸性季風氣候，夏季高溫多雨，冬季寒冷干燥，降水季節分配不均。地形的復雜多變造成了北京地區天氣狀況的多樣性。暴雨和強對流天氣是北京地區夏季(6-8月)主要的災害性天氣[8,9]，導致北京地區旱澇災害頻繁。如圖1所示為北京市1951-2015近65 a的降水過程。其中，年最大降水量為1 407 mm，發生在1957年；年最小降水量為267 mm，出現在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年月降雨資料，是經過該機構整編審查后發布的，故資料的可靠性基本可以得到保證。對資料的代表性進行分析，如圖2所示，繪制北京市年降雨量模比系數差積曲線，可以看出隨時間的增長，曲線變幅越來越小，年降雨系列在50年以上時，模比系數差積曲線穩定趨近于“1”，說明此系列資料具有一定的穩定性。

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

2.2　降水變化的多尺度分析

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

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

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

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

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

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

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

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

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

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

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

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

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

2.4　降水變化的趨勢分析

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

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

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

3　結　論

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

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

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

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

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

□

參考文獻：

[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]．物理學報，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]王金良，李宗軍．可用于氣候數據分析的ESMD方法[J]．氣候變化研究快報，2014，(3):1-5．

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

[8]北京志編委．北京志地質礦產水利氣象卷水利志[M]．北京：北京出版社，2000．

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

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

[11]魏鳳英．現代氣候統計診斷與預測技術[M]．2版.北京：氣象出版社，2007:36-104．

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

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

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

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

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