婁德君 潘昕濃 王冀 張雪梅 高振鐸

(1.齊齊哈爾市氣象局，黑龍江 齊齊哈爾 161006； 2.北京市氣象服務(wù)中心，北京 100089；3.北京市氣候中心，北京 100089； 4.哈爾濱市氣象臺，黑龍江 哈爾濱 150028 )

引言

對觀測事實的研究表明，氣候系統(tǒng)變化具有非平穩(wěn)特征[1-3]，而外部強迫隨時間的變化是導(dǎo)致非平穩(wěn)行為產(chǎn)生的根本原因[4-5]。慢特征是指信號源中隱含的變化最緩慢的成分，能夠表征信號源的固有屬性[6]。研究表明，慢特征分析方法(Slow Feature Analysis，SFA) 是從已知的非平穩(wěn)時間序列中提取緩變信息的有效方法，最近十余年，該方法在氣象領(lǐng)域得到了初步成功的嘗試[7-8]。根據(jù)Packard等[9]的狀態(tài)空間重構(gòu)理論以及Takens[10]的嵌入定理，王革麗等[11]在構(gòu)建預(yù)測方程時考慮了外強迫信息，建立了外強迫因子參與在內(nèi)的非平穩(wěn)時間序列預(yù)測模型。試驗結(jié)果表明，增加外強迫信息能夠有效地改善模型預(yù)測精度。張彬[12]、陳瀟瀟等[13]分別對500 hPa月平均高度場和小時氣溶膠濃度序列等進行了外強迫信息的提取，并將外強迫信息參與到預(yù)測建模中，結(jié)果均表明引入外強迫因子可以提高模型的預(yù)測技巧。

氣溫預(yù)測建模[14-15]時需考慮多種因素，而各因子之間存在著不同程度的相關(guān)，預(yù)測結(jié)果也會受到影響。以往研究表明，SFA方法對氣候系統(tǒng)的慢變信號提取避免了常規(guī)方法主觀選擇因子的不足，客觀性較強。因此，本文擬采用SFA方法對哈爾濱市氣溫序列進行研究，提取其慢特征信號，并構(gòu)建包含慢特征信號的預(yù)測模型，旨在增進對黑龍江省氣溫變化規(guī)律的了解，并為黑龍江省氣候預(yù)測提供一個新的思路。

1 資料與方法

1.1 資料來源

1961—2020年黑龍江省62個地面氣象站逐月氣溫資料來自黑龍江省氣候中心。相關(guān)分析表明黑龍江省氣溫具有較好的空間一致性，全年12個月各站點之間氣溫相關(guān)均較好，除個別站點在7月、8月、12月、5月與其余站相關(guān)略差外，大部分臺站與其余站點氣溫的相關(guān)都通過了95%置信度檢驗，相關(guān)系數(shù)的大值區(qū)主要位于黑龍江省中部地區(qū)(圖略)。哈爾濱市作為省會城市，且與全省其余61個臺站氣溫的相關(guān)系數(shù)均通過了99%以上信度檢驗，具有代表性。因此，本文選取哈爾濱市為黑龍江省代表站，進行氣溫的慢特征信號提取及預(yù)測研究。文中氣溫均值取1981—2010年30 a平均。

1.2 研究方法

主要方法有慢特征分析(Slow Feature Analysis，SFA)、小波變換、相關(guān)分析等。

1.2.1 SFA方法簡介

SFA方法是一種基于特征向量方法的無監(jiān)督算法，基本思想是對向量矩陣進行特征值分解，最終得到函數(shù)最優(yōu)解的集合[5,16]。慢特征分析方法步驟已有很多工作進行介紹[6,17-19]：

首先給定一個由非平穩(wěn)時間序列重構(gòu)的相空間X(t)，對其進行標(biāo)準(zhǔn)化處理

x(t)=[x1(t),x2(t),…,xn(t)]T

(1)

設(shè)置嵌入維數(shù)為m，時滯參數(shù)τ=1，對x(t)進行非線性擴展

(2)

(3)

構(gòu)建導(dǎo)函數(shù)空間，建立導(dǎo)函數(shù)的協(xié)方差矩陣并進行主成分分析，計算結(jié)果中選取最小特征值λj，對λj對應(yīng)的特征向量進行標(biāo)準(zhǔn)化即得到權(quán)重向量w1,…,wj，輸入函數(shù)變?yōu)間j(t)=wjz(t)。積分求得輸出信號

yj(t)=rwjz(t)+c

(4)

其中r與c均為常數(shù)。選取最小的特征值時，積分后得到變化最慢的信號分量，即慢特征信號。

1.2.2 SFA方法提取非平穩(wěn)系統(tǒng)中慢特征信號的能力

一維Logistic映射是一個非常簡單的混沌映射，Logistic映射的動力學(xué)方程為

xn+1=μxn(1-xn)[20]

(5)

其中μ為控制參數(shù)，μ∈(0,4]，當(dāng)μ∈[3.569945627,4]時，xn處于混沌區(qū)。

最大Lyapunov指數(shù)大于零通常作為系統(tǒng)存在混沌的判據(jù)。設(shè)定試驗的控制參數(shù)μ∈[3,4]，步長為0.001，將μ代入到映射方程中迭代10000步(圖略)可以看到，在μ取值區(qū)間內(nèi)，存在最大Lyapunov指數(shù)大于零，因此系統(tǒng)處于混沌狀態(tài)。

構(gòu)造帶有兩個時變參數(shù)at、bt的Logistic映射

xt+1=atxt(bt-xt)t=1,2,…,N

(6)

令參數(shù)遵從下式

at=3.5+0.3cos (4πt/500)

(7)

bt=1.02+0.01cos(15πt/500)

(8)

at和bt分別代表慢變和快變信號，atbt在3.236—3.914范圍內(nèi)變化，表明序列在周期區(qū)和混沌區(qū)變化，即該序列是非平穩(wěn)信號。

令初值x1=0.2，對式(6)迭代15000次，提取最后2000個數(shù)據(jù)作為試驗時間序列{xt}，令嵌入維數(shù)m=11，時滯參數(shù)τ=1，對{xt} 進行慢特征分析，提取SFA信號{yt}，各信號軌線變化見圖1。

圖1 Logistic映射的原始外強迫信號{at}(a)、原始外強迫信號{bt}(b)、非平穩(wěn)序列{xt}(c)和SFA信號{yt}(d)

{xt}為非平穩(wěn)時間序列，SFA方法提取到的慢特征信號{yt}與兩個原始真實外強迫信號均不相同，但其變化周期與慢變信號{at}的周期很相似， {yt}軌線上峰谷值附近呈鋸齒狀，可能是疊加了更小尺度的變化。采用Morlet小波變換方法[21]分析{yt}的時間平均功率譜(圖2a)，由圖可見{yt}功率譜有兩個峰值，均通過95%置信度檢驗(圖2a中虛線)。提取這兩個峰值點所對應(yīng)的尺度信號分量P1、P2，將標(biāo)準(zhǔn)化處理后的分量與兩個原始外強迫信號{at}、{bt}比較，可以看到平移了約半個周期長度的P1分量與{at}擬合較好(圖2b)，除兩端受邊界效應(yīng)影響曲線振幅差異較大外，兩曲線的中間時段擬合度較好，二者相關(guān)系數(shù)達到0.99。P2分量平移后與{bt}擬合度也較好(圖2c)，二者相關(guān)系數(shù)為0.98。

圖2 Logistic映射SFA信號的小波功率譜(a)、SFA信號慢變分量P1與原始外強迫信號{at}(b)、快變分量P2與原始外強迫信號{bt}(c)的比較

對Logistic映射構(gòu)造的非平穩(wěn)系統(tǒng)試驗表明，SFA方法可以從序列中提取出與原始外強迫信號擬合較好的慢特征信號。利用小波分析方法可以分離出SFA信號中的顯著周期分量，慢變和快變分量與原始外強迫信號都僅有位移和振幅差別，相關(guān)系數(shù)均超過0.98。

1.2.3 基于SFA方法的預(yù)測模型建立方法

王革麗等[11]在提取非平穩(wěn)系統(tǒng)慢特征信號的基礎(chǔ)上，結(jié)合狀態(tài)空間重構(gòu)理論、嵌入定理等，建立了非平穩(wěn)氣候預(yù)測模型：

{xi}、{Fi}為兩個非平穩(wěn)序列，{xi}是系統(tǒng)狀態(tài)隨時間的變化，{Fi}是{xi}的慢變信號。借助時滯參數(shù)τ將它們嵌入到m1+m2維的狀態(tài)空間中，得到狀態(tài)軌線

{xi,xi-k,…,xi-(m1-1)k;Fi,Fi-τ,…,Fi-(m2-1)k}i=1,2,…,N

(9)

式(9)中，m1為{xi}、m2為{Fi}的嵌入維數(shù)；N=n-(max(m1,m2)-1)k為狀態(tài)軌線上的點數(shù)。建立預(yù)測模型如下

(10)

問題即轉(zhuǎn)化為求解多元回歸函數(shù)問題。

2 結(jié)果分析

利用1961—2016年共672個月哈爾濱市月氣溫序列進行前期試驗，包括參數(shù)敏感性試驗、SFA信號分析等，對2017—2020 年48個月氣溫作獨立樣本預(yù)測。

2.1 對參數(shù)進行敏感性試驗

通常計算嵌入維數(shù)的方法是以嵌入定理為依據(jù)的，而氣候時間序列為非平穩(wěn)系統(tǒng)，該方法不再適用，因此需要依靠經(jīng)驗和試驗來分析SFA慢特征信號對參數(shù)選取的敏感性，最終確定適合研究的參數(shù)。選取哈爾濱市氣溫序列為實際非平穩(wěn)時間序列對m進行敏感性試驗，設(shè)置時滯參數(shù)τ=1，變化m的取值范圍，分別提取氣溫序列的慢特征信號。圖3給出了部分試驗結(jié)果，由圖可見，當(dāng)m取值為5、9、13時，提取的慢特征信號的變化趨勢和月際間變化都比較接近，只是在峰值和谷值附近存在微小差異。隨m值的增大，慢特征信號的振幅減小，曲線噪音減小，平滑度逐漸增大?？梢姡谝欢ǚ秶鷥?nèi)選取m值，提取出的慢特征信號差別不大。參考前人研究并結(jié)合敏感性試驗結(jié)果，本文對氣溫序列提取慢特征信號時設(shè)定嵌入維數(shù)m=11，時滯參數(shù)τ=1。

圖3 不同嵌入維數(shù)m=5(a)、m=9 (b)、m=13(c)、m=22 (d)、m=31(e)提取的哈爾濱市月氣溫慢特征信號比較

2.2 哈爾濱市氣溫的慢特征信號

依據(jù)敏感性試驗結(jié)果將參數(shù)設(shè)置為m=11，τ=1。利用SFA方法提取1961年1月至2016年12月哈爾濱市月平均氣溫距平序列的慢特征信號(圖4)。由圖4可見，慢特征分析方法提取的SFA信號更加平滑，但能較好地反映原序列的時間變化趨勢、高低轉(zhuǎn)換、極值等信息，進一步計算可知，提取的慢特征信號變率[22]對氣溫序列變率的貢獻率為57%。二者相關(guān)系數(shù)R=0.6，通過了99.9%信度的顯著性檢驗。

圖4 1961—2016年哈爾濱市月氣溫距平序列(a)及提取的SFA信號(b)變化

2.3 哈爾濱市氣溫SFA信號的時間尺度特征

對提取的哈爾濱市氣溫SFA信號進行Morlet小波變換，由圖5中SFA信號的小波功率譜可以看到，哈爾濱氣溫SFA信號具有4 a和7 a左右兩個較顯著的年際尺度周期振蕩，且2.0—7.8 a的周期均通過95%置信度的顯著性檢驗。與以往研究得到ENSO事件的2—7 a[23]年際尺度周期類似。相關(guān)分析表明，哈爾濱市氣溫SFA信號與前期3—12個月的NINO3.4區(qū)海溫指數(shù)顯著正相關(guān)，相關(guān)系數(shù)通過了95%信度檢驗。侯依玲等[24]研究也得到中國東北地區(qū)年平均氣溫一致型的時間系數(shù)存在顯著的2—7 a周期，指出中國東北北部地區(qū)氣溫受ENSO等大尺度氣候背景影響顯著。因此，哈爾濱氣溫的變化可能與前期赤道東太平洋海溫變化相關(guān)，但相應(yīng)的物理機制還需進一步探討。

圖5 1961—2016年哈爾濱市氣溫序列SFA信號的小波全局功率譜

2.4 哈爾濱市氣溫的預(yù)測試驗

在提取了哈爾濱市氣溫慢特征信息的基礎(chǔ)上，構(gòu)建包含慢特征信號的預(yù)測模型，并分析擬合和預(yù)測結(jié)果。在評價SFA信號預(yù)測能力的基礎(chǔ)上滾動預(yù)測近4 a哈爾濱市逐月氣溫距平。

試驗使用1961—2016年共672個月氣溫數(shù)據(jù)來進行建模及預(yù)測試驗，使用相同長度的150個子數(shù)據(jù)序列的平均預(yù)測結(jié)果來評價預(yù)測方程，主要評價指標(biāo)為方程擬合及預(yù)測的相關(guān)系數(shù)和均方根誤差(RMSE)。

每個子序列前面40 a的數(shù)據(jù)用于提取SFA信號并建立預(yù)測函數(shù)，預(yù)測時間步長為1—12個月。時滯參數(shù)τ取為1，{xt}為氣溫距平序列，嵌入維數(shù)m1取為3—8，{αt}為慢特征信號，嵌入維數(shù)m2取為0—3。m2=0代表預(yù)測方程中不包含慢特征信號，即平穩(wěn)性模式。m2=1—3代表預(yù)測方程中嵌入1—3個維數(shù)的慢特征信號，簡稱慢特征模式。

由圖6a可以看到，m2=0平穩(wěn)性模式擬合的相關(guān)系數(shù)在預(yù)測時間步長為1—2，即超前1—2個月預(yù)測時，通過了95%置信度的檢驗。時間步長大于2時，相關(guān)系數(shù)保持在0.2左右，沒有通過95%置信度檢驗。當(dāng)m2=1—3，即加入了慢特征信號后，前1—5步擬合相關(guān)系數(shù)明顯增大，且隨嵌入維數(shù)m2增大相關(guān)系數(shù)略增大。除了m2=1在第4—5步時相關(guān)系數(shù)值在0.4—0.6之間外，其余情況第1—5步的相關(guān)系數(shù)值都超過了0.6。而從第6步開始，即預(yù)測超前6個月以上時，擬合相關(guān)系數(shù)急劇下降，接近95%置信度臨界線。在全部時間步長里，慢特征模式的相關(guān)系數(shù)都高于平穩(wěn)性模式，步長為1—5時，慢特征模式相關(guān)系數(shù)偏高明顯。由圖6b可見，預(yù)測的相關(guān)系數(shù)m2=0時在全部時間步長均未通過95%置信度檢驗。m2=1—3時曲線形態(tài)在預(yù)測時間步長為1—5步時和擬合曲線接近，相關(guān)系數(shù)都通過95%信度檢驗，但值略低于擬合結(jié)果，第6步開始相關(guān)系數(shù)值急劇下降至0值附近，明顯低于圖6a的0.25，即第6步以后可預(yù)報性較小。在圖6a和圖6b中，相關(guān)系數(shù)最大值均出現(xiàn)在預(yù)測時間步長為第5步時。

圖6 氣溫預(yù)測模型擬合的(a)、預(yù)測的(b)相關(guān)系數(shù)和擬合的(c)和預(yù)測的(d)均方根誤差

擬合的均方根誤差(圖6c)在m2=0時明顯偏大，誤差值在1.7—1.9之間變化，預(yù)測步長為1時，誤差最小，第2—3步明顯上升，之后穩(wěn)定維持。m2=1—3時，即加入慢特征信號后，在預(yù)測步長為1—5步時，均方根誤差明顯偏小，誤差值在1.4—1.5之間變化，第5步時誤差有明顯下降，之后誤差明顯增加，第7步以后接近常值，但略小于m2=0的誤差。隨嵌入維數(shù)m2增加，均方根誤差略有減小。預(yù)測的均方根誤差(圖6d)形態(tài)與擬合的均方根誤差相似，但誤差值整體均略大于擬合的誤差值。擬合和預(yù)報結(jié)果均在時間步長為第5步時，均方根誤差最小。

試驗結(jié)果表明，引入慢特征信號后建立的哈爾濱市月氣溫預(yù)測模型在預(yù)測時間步長為1—5，即提前1—5個月預(yù)測時，相關(guān)系數(shù)明顯高于平穩(wěn)性模型，預(yù)測誤差也明顯小于平穩(wěn)性模型。這表明在一定的預(yù)測時間步長內(nèi)，引入慢特征信號的預(yù)測模型可以提高預(yù)測能力和精度。嵌入維數(shù)m2=3，預(yù)測時間步長為5時，預(yù)測效果最優(yōu)。

2.5 哈爾濱市月氣溫預(yù)測

依據(jù)2.4試驗方案，對哈爾濱市2017年1月至2020年12月共48個月氣溫進行滾動預(yù)測。設(shè)定相空間嵌入維數(shù)m1=8，m2=3，時間步長k=5。逐月氣溫預(yù)測及與實況對比結(jié)果見圖7。預(yù)測結(jié)果統(tǒng)計表明，平穩(wěn)性模型m2=0的預(yù)測同號率為0.40左右(19/48)，預(yù)測平均均方根誤差為1.7 ℃，而慢特征模式m2=3預(yù)測同號率為0.71(34/48)，預(yù)測平均均方根誤差為1.5 ℃?？梢?，引入慢特征信號的預(yù)測模型效果優(yōu)于平穩(wěn)性模型。

圖7 氣溫獨立樣本預(yù)測與實況的比較

3 結(jié)論與討論

(1)慢特征分析方法可以從Logistic映射時間序列中提取出與原始外強迫信號擬合較好的慢特征信號，其慢變和快變分量與原始外強迫信號僅有位移和振幅差別。

(2)SFA方法可以有效地提取哈爾濱市氣溫序列中的慢特征信號。提取的SFA信號能夠反映原氣溫序列曲線的變化趨勢、極值等信息。

(3)擬合和預(yù)測試驗表明，引入SFA信號后的哈爾濱市氣溫預(yù)測模型可以在一定程度上提高預(yù)測能力和精度。48個月獨立樣本滾動預(yù)測試驗也表明引入慢特征信號后預(yù)測效果優(yōu)于平穩(wěn)性模型。

(4)慢特征分析方法中參數(shù)的選擇，是通過試驗嘗試并參考以往文獻主觀選取的，缺乏客觀性判斷指標(biāo)，如何更加合理地選取嵌入維數(shù)等參數(shù)還需要進一步研究。利用SFA方法提取慢特征信號只是一種數(shù)學(xué)算法，并不能解釋提取的慢特征信號所代表的物理含義。功率譜分析能夠得到慢特征信號的顯著周期尺度，但其所代表的外部影響因子，也需要相應(yīng)的物理機制來佐證。黑龍江省地處中高緯，氣溫變化的成因非常復(fù)雜，除受外強迫因子影響外，大氣環(huán)流因子的影響也很重要。提取的慢特征信號與環(huán)流因子之間的作用關(guān)系，二者對氣溫變率的貢獻等還需要結(jié)合更多的診斷分析方法來探討。