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

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

引言

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

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

1 資料與方法

1.1 資料來(lái)源

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

1.2 研究方法

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

1.2.1 SFA方法簡(jiǎn)介

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

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

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

(1)

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

(2)

(3)

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

yj(t)=rwjz(t)+c

(4)

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

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

一維Logistic映射是一個(gè)非常簡(jiǎn)單的混沌映射，Logistic映射的動(dòng)力學(xué)方程為

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

(5)

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

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

構(gòu)造帶有兩個(gè)時(shí)變參數(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分別代表慢變和快變信號(hào)，atbt在3.236—3.914范圍內(nèi)變化，表明序列在周期區(qū)和混沌區(qū)變化，即該序列是非平穩(wěn)信號(hào)。

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

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

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

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

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

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

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

{xi}、{Fi}為兩個(gè)非平穩(wěn)序列，{xi}是系統(tǒng)狀態(tài)隨時(shí)間的變化，{Fi}是{xi}的慢變信號(hào)。借助時(shí)滯參數(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)軌線上的點(diǎn)數(shù)。建立預(yù)測(cè)模型如下

(10)

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

2 結(jié)果分析

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

2.1 對(duì)參數(shù)進(jìn)行敏感性試驗(yàn)

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

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

2.2 哈爾濱市氣溫的慢特征信號(hào)

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

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

2.3 哈爾濱市氣溫SFA信號(hào)的時(shí)間尺度特征

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

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

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

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

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

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

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

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

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

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

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

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

圖7 氣溫獨(dú)立樣本預(yù)測(cè)與實(shí)況的比較

3 結(jié)論與討論

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

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

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

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