奇異譜分析奇異譜分析最大熵預(yù)報方法

奇異譜分析奇異譜分析最大熵預(yù)報方法

ssa—引言異常光譜分析的最大熵預(yù)處理方法（ssa-mem）是最早用于氣象研究的東西海岸和格里爾。目前,California大學的AmiraSaunders、MichaehGhil等用SSA—MEM做～Νino3區(qū)海溫距平(SSTA)和南方濤動指數(shù)(SOI)預(yù)報。SSTA本身的功能較多,MEM也是時間序列建模和譜估計的一個重要方法,但要在有限篇幅內(nèi)完整地介紹它們很困難。因此,詳細了解SSA可參考文獻,了解MEM可參考一般的時間序列分析與應(yīng)用書籍。這里側(cè)重介紹業(yè)務(wù)人員關(guān)心的應(yīng)用SSA-MEM建立預(yù)報方法的步驟。1原序列及一般序列中時間主成分的分析SSA—MEM是對一維時間序列做延伸預(yù)報的方法。要對一個序列做預(yù)報,首先必須認識該序列變化的主要規(guī)律。SSA是分析包含噪音的有限長序列的有效方法,它的主要功能是降低序列中的噪音、識別傾向成分和周期成分、與最大熵譜相結(jié)合估計序列的功率譜、消除高頻段的虛假譜峰。這里僅介紹與SSA—MEM有關(guān)的功能。SSA的分析對象是一維時間序列x1,x2,x3,…xn,要求中心化,平均值為零,也就是距平或標準化距平。我們的目的是要了解它的時間演變結(jié)構(gòu),因此,把它在時間上滯后排列X=[x1x2?xi+1?xΝ-Μ+1x2x3?xi+2?xΝ-Μ+2????xΜxΜ+1?xi+Μ?xΝ]≡[X1?0X1?1?X1?Ν-ΜX2?0X2?1?X2?Ν-Μ???XΜ?0XΜ?1?XΜ?Ν-Μ](1)X的第i個狀態(tài)用向量表示為Xi=[xi+1xi+2?xi+Μ]≡[X1?iX2?i?XΜ?i]i=0,1,...?Ν-Μ(2)共N-M+1個狀態(tài),X稱相空間中的軌跡矩陣,矩陣X中的元素與原一維時間序列對應(yīng)關(guān)系為:Xji=xj+i(3)把(1)式看作M個變量N-M+1個樣本。(1)式的變量相互間的協(xié)方差是原序列xi不同滯后的自協(xié)方差。對其做主成分分析,步驟為先計算滯后自協(xié)方差矩陣TxΤx=[C(0)C(1)C(2)?C(Μ-2)C(Μ-1)C(1)C(0)C(1)?C(Μ-3)C(Μ-2)C(2)C(1)C(0)?C(Μ-4)C(Μ-3)?????C(Μ-2)C(Μ-3)C(Μ-4)?C(0)C(1)C(Μ-1)C(Μ-2)C(Μ-3)?C(1)C(0)](3)Tx是Toeplitz矩陣,也是實對稱矩陣。它的元素C(j)用x1,x2,…xN序列估計,即C(j)=1Ν-jΝ-j∑t=1xtxt+jj=0,1,2,??Μ-1(4)然后用雅可比子程序計算Tx的特征向量Ek和特征值λk,滿足的方程為TxEk=λkEk,k=1,2,…,M(5)其中Ek是Tx與λk對應(yīng)的歸一化特征向量,λk從大到小排列,λ1>λ2>…>λM≥0。再計算X的狀態(tài)向量Xi在Ek上的投影aki=Xi?Ek=Μ∑j=1XjiEkj=Μ∑j=1xi+jEkj0≤i≤Ν-Μ(6)把以上運算與空間場的EOF分解相對照,這里的特征向量Ek反映x序列中的時間演變型,稱為時間EOF(簡稱T-EOF)。aki是Ek所反映的時間演變型在原序列的xi+1,xi+2,…,xi+M時段的權(quán)重,稱為時間主成分(簡稱T-PC)。得到T-EOF和T-PC以后,可以與傳統(tǒng)的空間EOF分析一樣,由前幾個重要特征向量和時間系數(shù)重建原來的分析對象,重建成分簡稱RC(ReconstructionComponents)。由第k個T-EOF和T-PC重建xi的成分記為xki,公式為xki={1ΜΜ∑j=1aki-jEkjΜ≤i≤Ν-Μ+11ii∑j=1aki-jEkj1≤i≤Μ-11Ν-i+1Μ∑j=i-Ν+Μaki-jEkjΝ-Μ+2≤i≤Ν(7)所有RC之和等于原序列xi=Μ∑k=1xkii=1,2??Ν(8)截取前k個貢獻大的成分近似表示原序列?xi=k∑k=1xkii=1,2??Ν(9)?xi序列一般集中了原序列xi中時間尺度大的成分,不太嚴格地說,也是相對容易做預(yù)報的成分。k取多大,可通過反復(fù)做運算比較確定,也可參照下面的分析確定。SSA識別某個RC是否屬于傾向成分采用Kendall非參數(shù)檢驗。考察第k個RC時,計數(shù)滿足xki<xkj的指標(i,j)數(shù)Kr,統(tǒng)計量τ=4ΚrΝ(Ν-1)-1(10)當xk不是傾向成分的原假設(shè)成立時,τ服從均值為零,均方差為S=√2(2Ν+5)9Ν(Ν-1)(11)的正態(tài)分布。因此,若取信度α=0.05,當樣本τ值落在-1.96S和1.96S以外時,拒絕原假設(shè),認為第k個RC是傾向成分。τ>1.96S和τ<-1.96S分別對應(yīng)上升和下降傾向。對一些實際氣象序列和人造序列計算表明,這個檢驗?zāi)茏R別出傾向成分,但標準不高,往往會識別出幾個傾向成分。這時應(yīng)注意再從中挑選τ絕對值最大的確定1～2個RC為傾向成分,而滿足以上條件的其他RC可能是長周期成分。往往由于原序列xi不夠長,長周期成分在時間t=1～N內(nèi)不是整數(shù)個周期循環(huán),在兩端不完整的循環(huán)部分呈上升或下降趨勢,就可能使τ值達上述標準而被識別成傾向成分。根據(jù)SSA原理,當原序列中存在一個周期成分時,SSA將得到一對RC,它們的特征值接近相等(因此一定是相鄰序號的一對,不失一般性可設(shè)這對特征成分的充號為k和k+1),對應(yīng)的一對T-EOF和T-PC分別正交。原序列中的一個周期成分是滿足這些條件的一對RC之和?，F(xiàn)舉例說明,設(shè)有一個周期成分xt=Acos(ωt+φ0)(12)SSA相當于從(1)式矩陣出發(fā)做主成分分析,其中的Xji與原序列的關(guān)系為Xji=xj+i,記時間取樣間隔為△t,則Xji=Acos(ω(j+i)△t+φ0)=Acos(ωi△t+φ0)cos(ωj△t)-Asin(ωi△t+φ0)sin(ωj△t)(13)與主成分的展開式Xji=∑kakiEkj(14)相對照,滯后空間(j為指標)的函數(shù)是cos(ωj△t)和sin(ωj△t),它們的頻率相等,若將其歸一化,即是一對正交的T-EOF(Ek和Ek+1),樣本序號空間(i為指標)的函數(shù)是Acos(ωi△t+φ0)和-Asin(ωi△t+φ0),是一對時間主成分,它們也正交。由于特征值等于主成分的方差,而Acos(ωi△t+φ0)和-Asin(ωi△t+φ0)的方差都等于A2/2,故λk=λk+1。然而,上述例子是理想化的情況,對離散的有限長序列做SSA,即使是周期信號也不能得到嚴格相等的λk與λk+1,T-EOF間和T-PC間雖總是正交的,但T-EOFx和T-EOFk+1有嚴格相等的頻率也不現(xiàn)實。為此,Vautard和Ghil提出了補充和修正的判據(jù):①λk和λK+1接近相等。②T-EOFk和T-EOFk+1有相近的頻率,判斷方法是對T-EOFk和T-EOFk+1做付里葉變換,Ek的付氏變換是～Ek(f)=Μ∑j=1Ekjei2πfi(15)式中i=√-1,因此～Ek(f)是復(fù)數(shù)。在f=0.0～0.5間的500個等距f上計算付氏變換的模方|～Ek(f)|2和|～Ek+1(f)|2,找出|～Ek(f)|2和|～Ek+1(f)|2達最大值的頻率,分別記為fk和fk+1。記δfk=|fk-fk+1|,δfk應(yīng)很小。因為對純紅噪音過程有關(guān)系δfk=1/(2M),故至少要求δfk<1/(2M)。由于我們想識別周期振蕩成分,故尚需加強此標準,可要求2Mδfk<0.75。③|～Ek(f)2|和|～Ek+1(f)|2足夠大。若原序列中頻率介于fk和fk+1之間的頻率f*的周期振蕩完全由這一對RC表示出來,則可導出1Μ[|～Ek(f*)|2+|～Ek+1(f*)|2]=1(16)在實際情況下,一對成分不可能與一孤立的頻率完全對應(yīng),所以取1Μ[|～Ek(f*)|2+|～Ek+1(f*)|2]>23為標準,這意味著原序列xi中頻率為f*振蕩的方差至少有2/3被這對重建的成分表示出來。滿足上述條件的一對RC之和是原序列中的一個周期振蕩成分。文獻對SOI和～Νino3區(qū)的SSTA做SSA,得出第1、2RC成對,周期是準2年;第3、4RC成對,周期是準4年。2譜估計軟件“推動自回歸”xt序列的最大熵譜是一個p階自回歸模型的功率譜,?Gx(f)=Ρp|1+p∑k=1ake-i2πfk|2(17)這里不介紹最大熵譜的概念和算法,只說明SSA-MEM中利用到的部分,并且使用文獻中的公式和程序。(17)式中的ak是以下自回歸模型的回歸系數(shù)xt+a1xt-1+a2xt-2+…+apxt-p=εt(18)最優(yōu)階數(shù)p的確定可用最終預(yù)報誤差(FPE)準則或赤池準則(AIC)。作為預(yù)報方程,寫為?xt=φ1xt-1+φ2xt-2+?+φpxt-p(19)顯然φk=-ak,φk可以利用文獻所附的“極大墑譜估計的伯格算法”程序(22MEBURG)得到,即利用該程序計算?Gx(f)的ak。使用該程序應(yīng)注意兩點。首先它需讀入資料序列及一些參數(shù),各參數(shù)的意義從程序的說明中容易看出。其中讀入一個<40的最高階數(shù)(變量名MLAST)。程序的計算結(jié)果是這個給定的最高階自回歸模型的功率譜,不是最優(yōu)階數(shù)的。相應(yīng)存放自回歸系數(shù)的數(shù)組最后結(jié)果也是最高階自回歸模型的系數(shù),因此需變更一下。最簡單的做法是,第一次給最高階數(shù)38或39,執(zhí)行一次,除得到最高階模型的功率譜外,它還輸出最優(yōu)自回歸的階數(shù)(變量名NAPOPT),然后把所得到的最優(yōu)階數(shù)作為最高階數(shù)讀入,再執(zhí)行一次程序。當然,把該程序修改成一次執(zhí)行即可得到最優(yōu)階自回歸模型的功率譜及回歸系數(shù)也不困難。其次,對照公式(17)閱讀該程序可以看出,程序中G數(shù)組的元素與ak的關(guān)系是ak=G(k+1),因此,預(yù)報方程(19)中的系數(shù)是φk=-G(k+1)。G(k+1)是G數(shù)組的第k+1個元素)。這樣,對時間序列可以建立一個最優(yōu)階數(shù)的自回歸預(yù)報模型(AR)。在SSA-MEM中,AR預(yù)報模型不是直接用于x序列,而是用于時間主成分aki序列。優(yōu)點在于,根據(jù)(6)式,aki是原序列x的加權(quán)滑動平均,加權(quán)系數(shù)是Ek的分量,所以aki(0≤i≤N-M)是原x序列的一種濾波結(jié)果,較之原序列x和aki序列變化單調(diào)些,具有更多的可預(yù)報性。對前幾個aki序列(k=1,2,…K)分別建立各自的AR模型。可把aki向后預(yù)報τ步,得?akΝ-Μ+1、…、?akΝ-Μ+τ,再把預(yù)報出的時間主成分與SSA得到的T-EOF相結(jié)合即得原序列的預(yù)報?xN+1、?xN+2、…、?xN+τ。計算公式類似于(7)的第3式,或參照空間EOF分析中由特征向量和時間系數(shù)恢復(fù)原變量的公式。3氣候和持續(xù)性預(yù)報標準還算例2SSA-MEM用于SOI和～Νino3區(qū)SSTA預(yù)報的效果還可以。是否可用于其他氣候要素預(yù)報應(yīng)先實驗,并與氣候預(yù)報和持續(xù)性預(yù)報比較。氣候預(yù)報是用氣候平均值作為預(yù)報值,因為x序列是距平,所以就是用零作為該序列的預(yù)報值。持續(xù)性預(yù)報是用xN作為以后的預(yù)報?xN+1,?xN+2,?xN+3,…。以下舉例說明效果估計方法。假如已有x序列的歷史資料長度為606個月,則可取N=561,562,…,600,共做40次實驗,每次預(yù)報以后6個月的。這樣,SSA-MEM預(yù)報、氣候預(yù)報、持續(xù)性預(yù)報都可以與實況相比較,計算誤差平方和。ef(τ)g=(?xΝg+τ-xΝg+τ)2ec(τ)g=(0.0-xΝg+τ)2ep(τ)g=(xΝg-xΝg+τ)2上式分別是3種預(yù)報的誤差平方,g表示第幾次實驗,g=1～40,Ng是第g次試驗的序列長度,Ng=561,562,…,600,然后求均方根誤差。RΜEf(τ)=√14040∑g=1ef(τ)g