基于Lorenz模型的集合預報與單一預報的比較研究?

梁丁顧斌丁瑞強李建平鐘權加

1)(南京信息工程大學空間天氣研究所,南京 210044)

2)(中國科學院大氣物理研究所,大氣科學和地球流體力學數(shù)值模擬國家重點實驗室,北京 100029)

3)(南京信息工程大學物理與光電工程學院,南京 210044)

4)(北京師范大學全球變化與地球系統(tǒng)科學研究院,北京 100875)

1 引 言

大氣是一個復雜的非線性系統(tǒng),任何初始的微小誤差都會導致預報結果存在很大的不確定性[1?3].Epstein[4]和Leith[5]提出集合預報解決“單一確定性”預報存在“不確定性”的問題.其基本思想是:針對大氣初始狀態(tài)的不確定性,選取合適的擾動集合,疊加于預報初始分析場上,然后分別向前積分制作預報,得到預報集合.大量研究和業(yè)務應用表明,集合預報相對單一預報而言,能有效減小預報誤差,改善預報水平,提高預報技巧[6,7].

近年來,集合預報的研究重點主要是:1)初值不確定性;2)模式不確定性;3)集合預報產品釋用.在不考慮模式不確定性的理想情況下,集合預報成功的關鍵在于初始擾動的產生.圍繞初始擾動的生成,眾多的理論方法被提出,如蒙特卡羅法[5]、時間滯后平均法[8]、增長模繁殖(bred growing mode,BGM)法[9,10]、奇異向量(singular vector,SV)法[11?13]、觀測擾動法[14,15]、條件非線性最優(yōu)擾動法[16,17]、集合卡爾曼濾波(ensemble Kalman fi lter,EnKF)法[18]、集合變換卡爾曼濾波(ensemble transform Kalman fi lter,ETKF)法[19]以及新近發(fā)展的非線性局部Lyapunov向量(nonlinear local Lyapunov vector,NLLV)法等[20?22]. 以上述理論方法為基礎,國內外學者開展了許多研究工作.例如,Anderson[23]基于Lorenz63理想模型,評估了BGM,SV以及隨機擾動等方法的預報表現(xiàn),指出了隨機擾動方法比BGM和SV方法具有更高的預報技巧;封國林等[24]利用Lorenz63模型,對集合預報可能的物理基礎進行了探討;Bowler[25]基于Lorenz96模型的研究結果表明,EnKF方法的預報效果比BGM和SV方法更優(yōu).值得注意的是,這些工作主要圍繞初始擾動產生方法的比較,雖然取得了許多有意義的結果,但多數(shù)研究為大樣本平均的結果,較少涉及預報技巧對個例依賴性問題的研究.集合預報的預報表現(xiàn)對具體實驗個例的敏感性尚需進一步研究.此外,雖然人們對集合預報初始擾動的生成及其在業(yè)務預報中的應用等方面的研究較多[26?31],但對集合平均減小預報誤差的物理本質與基礎仍需進一步探討.

本文以Lorenz63和Lorenz96理想模型為例,利用NLLV和BGM方法,選取模型的不同狀態(tài)為實驗個例,研究了集合預報對不同實驗個例的預報表現(xiàn),將結果與單一預報進行對比,并試圖從概率分布演化的角度,初步分析集合平均具有較小預報誤差可能原因,希望能為集合預報理論和預報技巧的改善提供一定的依據(jù).

2 方法與物理模型

2.1 集合預報初始擾動集合產生

BGM方法是目前生成集合預報初始擾動集合的流行方法之一.該方法被美國國家環(huán)境預報中心(National Center for Environmental Prediction,NCEP)長期采用.其借鑒了分析循環(huán)中誤差增長過程,疊加任意隨機擾動,通過模式循環(huán)積分使其得到“繁殖”,當該擾動在某種意義上達到飽和時,擾動的大小和分布可被視為實際分析誤差快速增長模的估計.BGM方法的具體繁殖循環(huán)流程如圖1(a)所示:初始時刻,在分析場上疊加一組大小相等、方向隨機的擾動,并對受擾狀態(tài)獨立積分.在每個繁殖循環(huán)的末尾,將擾動的大小尺度化到初始時刻的擾動大小,并疊加到新的分析場上,進行下一個繁殖循環(huán).如此重復,直至預報時刻,就可以獲得一組BGM擾動.

圖1 (a)增長模繁殖擾動[10]、(b)非線性局部Lyapunov擾動生成[20]示意圖以及(c)集合預報流程圖Fig.1.Schematic of(a)bred growing modes[10],(b)nonlinear local Lyapunov[20]perturbations and(c)the flow chart of ensemble forecasting.

理論上BGM方法產生的擾動集合能很好地反映分析誤差中的增長成分,但在基本流的演化過程中,任何初始隨機擾動都會自然地向相空間中增長最快的主導線性局部Lyapunov向量(linear local Lyapunov vector,LLV)方向偏轉[9,10].因此,由BGM方法產生的擾動集合往往存在較高的相似性,從而影響其全局正交性[19,25],導致集合離散度不夠.針對上述BGM 方法的局限性,丁瑞強、李建平等[32?35]提出了在非線性理論框架下生成初始擾動的NLLV方法.該方法將非線性局部Lyapunov指數(shù)(nonlinear local Lyapunov exponent,NLLE)[20?22]拓展至多維,代表正交的擾動增長(縮小)方向.NLLV方法在保留BGM方法的優(yōu)點,即對參考和擾動狀態(tài)直接積分,獲取擾動集合的基礎上,加入正交化過程,保證擾動嚴格正交,因而能更好地抓住分析誤差的特征.NLLV的繁殖循環(huán)流程如圖1(b)所示:基本繁殖過程與BGM方法類似,不同的是在每個繁殖循環(huán)的末尾,對擾動進行正交操作,得到一組正交的擾動,再疊加到新的分析場上.如此重復,直至預報時刻.需要注意的是,在進行正交操作時,要先確定增長最快的擾動,記為主導非線性局部Lyapunov向量(leading nonlinear local Lyapunov vector,LNLLV),其余NLLV擾動可在LNLLV約束下依次獲得.Feng等[20,21]的研究結果表明,NLLV較BGM方法具有更高的預報技巧;與ETKF方法相比,兩者預報技巧接近,但NLLV計算所需資源要遠小于ETKF.本文采用BGM和NLLV方法生成集合預報初始擾動集合,進行相關預報實驗.

2.2 Lorenz模型與實驗設計

Lorenz63模型的控制方程如下:

其中,σ=10,γ=28.0,b=8/3,分別為Prandtl數(shù)、Rayleigh數(shù)以及表示與對流尺度相聯(lián)系的參數(shù).數(shù)值積分采用四階龍格-庫塔方法,積分步長取0.01,記作0.01 tus.在該非線性模型混沌吸引子上,每間隔0.05 tus選取一個狀態(tài),作為實驗個例的初始態(tài),共計選取N1(N1=104)個真實狀態(tài),記為Xt63.

Lorenz96模型作為大氣的低階近似,相較Lorenz63模型更適合于研究大氣,已經被廣泛應用于大氣預測和資料同化的研究中[36].Lorenz96模型對初值異常敏感,且非線性作用強,其動力框架可描述為

其中,Xi為模型變量,F為定常強迫.本文的實驗設計中,模型維數(shù)等于40,即i=1,2,………,40,強迫項F為8,數(shù)值積分采用四階龍格-庫塔格式,積分步長為0.05.對于Lorenz96模型而言,從誤差增長程度來看,間隔0.2時間單位大致對應1 d[37].因此本文Lorenz96模型的時間單位統(tǒng)一用天(d)表示.在該非線性模型混沌吸引子上每間隔0.25 d選取N2(N2=N1)個不同狀態(tài),記為Xt96,作為Lorenz96模型實驗個例的初始態(tài).

為了和實際預報情況接近,我們利用EnKF方法對Xt63,Xt96進行同化,生成對應的分析場Xa63,Xa96.作為比較,將從分析場Xa63,Xa96開始的預報視為單一預報.

Lorenz63模型任意實驗個例i(i=1,2,………,N1)的集合預報實驗設計如圖1(c)所示:在T1時刻的分析場上疊加M(M=2)個初始隨機擾動,擾動的模與真實分析誤差|Xa63i?Xt63i|一致,分布服從[?1,1]上均勻隨機分布;從疊加擾動后的分析場出發(fā),分別進行周期為0.1 tus,總時長為0.4 tus的繁殖循環(huán),到T0時刻,產生M個NLLV和BGM擾動;在T0預報時刻的分析上分別疊加和扣除M個擾動后得到2M個擾動場,將擾動場和分析場分別向前積分5 tus,對所有的預報成員采取等權平均,得到集合平均預報.

Lorenz96模型集合預報實驗設計與Lorenz63模型類似:在T1時刻的分析場上疊加M(M=5)個初始隨機擾動,擾動的模與真實分析誤差|Xa96i ? Xt96i |相等,分布服從[?1,1]上均勻隨機分布;從疊加擾動后的分析場出發(fā),分別進行周期為1 d,總時長為5 d的繁殖循環(huán),到T0時刻,產生M個NLLV和BGM擾動;在T0預報時刻的分析上分別疊加和扣除M個擾動后得到2M個擾動場,將擾動場和分析場分別向前積分12 d,對所有預報成員采取等權平均,得到集合平均預報.

3 結果與討論

3.1 集合平均預報的整體表現(xiàn)

首先以均方根誤差(root mean square error,RMSE)[38]和型異常相關(pattern anomaly correlation, PAC)[38]為標準,通過與單一預報的比較,評估Lorenz63模型和Lorenz96模型集合平均的整體預報表現(xiàn).圖2給出了所有實驗個例的平均RMSE和PAC演化過程.比較發(fā)現(xiàn),對Lorenz63和Loren96模型而言,隨著時間推移,NLLV和BGM兩種集合平均的預報技巧較單一預報都有明顯的改善,而NLLV比BGM表現(xiàn)更優(yōu).在預報初期(圖2((a),(b))中0.6 tus 前;(圖2(c),(d))中5 d前),NLLV和BGM集合平均與單一預報結果接近,隨著時間的推移,曲線開始出現(xiàn)分離,意味著集合平均的優(yōu)勢開始顯現(xiàn).在預報初期,誤差主要以線性增長為主,使得疊加的正負擾動預報對相互抵消.但隨著預報時間的延長,非線性作用增強,誤差逐漸進入非線性增長階段,集合平均的非線性濾波作用減小了預報誤差[10,39,40].

圖2 Lorenz63模型104個個例的平均RMSE(a)和平均PAC(b),Lorenz96模型104個個例的平均RMSE(c)和平均PAC(d)隨時間的變化(黑色實線,單一預報;藍色實線,BGM集合預報;紅色實線,NLLV集合預報)Fig.2.(a)Mean RMSE and(b)mean PAC of 104cases of the Lorenz63 model and(c)mean RMSE and(d)mean PAC of 104cases of the Lorenz96 model as a function of lead time(black solid curve,single forecast;blue solid curve,BGM ensemble forecast;red solid curve,NLLV ensemble forecast).

3.2 個例集合預報表現(xiàn)

3.1節(jié)的結果表明,集合平均的預報技巧高于單一預報.但對具體實驗個例而言,集合預報是否一定優(yōu)于單一預報仍不確定.為此,我們從預報誤差的角度,分析上述個例的集合預報表現(xiàn).

Lorenz63模型(圖3(a)—(b)),Lorenz96模型(圖3(d)—(f))所有實驗個例的NLLV集合平均和單一預報誤差如圖3所示.兩模型的結果基本一致,預報初期,集合平均和單一預報的預報誤差主要集中在“對角線”附近,意味著多數(shù)實驗個例,其集合平均和單一預報誤差大致相當.隨時間的增加,單一預報誤差大于集合平均的實驗個例數(shù)逐漸增多,多數(shù)實驗個例落入單一預報一側(BGM的結果與NLLV一致).

進一步集合平均誤差小于單一預報誤差的實驗個例占總實驗數(shù)的比例(k)如圖4所示.在預報開始階段,k≈50%.隨時間推移,k逐漸增大,集合平均的優(yōu)勢愈發(fā)顯著.Lorenz63模型第5 tus時刻的NLLV和BGM集合平均預報的k值分別為65%和62%.Lorenz96模型第12 d時刻的NLLV和BGM集合平均預報的k值分別為95%和86%.值得注意的是,對所有預報時刻,NLLV集合平均占優(yōu)的百分比都較高于BGM,這可能是由于NLLV的嚴格正交,導致其能捕捉更多的誤差分量,從而更準確地刻畫初始誤差向量[20,21].正如圖5所示,對于Lorenz63模型而言,當γ>24.74時,系統(tǒng)表現(xiàn)出不同復雜程度的混沌行為,計算表明,NLLV的這一優(yōu)勢不隨γ值的變化而變化.

圖3 Lorenz63和Lorenz96模型104個個例在第(a)1 tus,(b)3 tus,(c)5 tus和在第(d)1 d,(e)7 d,(f)12 d時刻的NLLV集合平均和單一預報的預報誤差Fig.3.Forecast error of NLLV ensemble mean and the single forecast of 104cases of Lorenz63 and Lorenz96 model at time(a)1 tus,(b)3 tus,(c)5 tus and at time(d)1 d,(e)7 d,(f)12 d,respectively.

圖4 Lorenz63模型(a)和lorenz96模型(b)集合平均預報比單一預報誤差小的實驗個例百分比隨時間的演變(點虛線,BGM集合平均;點實線,NLLV集合平均)Fig.4.Ratio of the cases of the NLLV(dot-dashed line)and BGM(dot-solid line)ensemble mean forecast with smaller error compared with the single forecast as a function of time from(a)Lorenz63 model and(b)Lorenz96 model.

圖5 Lorenz63模型NLLV(實線)和BGM(虛線)104個實驗個例在預報第5 tus時的k值隨γ的變化Fig.5.Ratio of the NLLV(solid line)and BGM(dashed line)ensemble mean forecast with smaller error compared with the single forecast at 5 tus as a function ofγof Lorenz63 model.

3.3 吸引子概率分布特征

吸引子的概率分布反映了系統(tǒng)的長期行為特征.吸引子上不變的概率分布是其基本特征[41].圖6給出了Lorenz96模型的變量X(X為任意變量)在不同長度的積分序列下的概率分布.很明顯可以看出,隨著演化時間趨于無窮,變量X的空間概率分布趨向不變,這意味著系統(tǒng)軌跡會以不變的概率在吸引子上演化.

圖7給出了不同預報時刻Lorenz96模型變量X的集合平均預報狀態(tài)和單一預報狀態(tài)的概率分布.對單一預報狀態(tài)而言,其概率分布在各預報時刻都與真實狀態(tài)的概率分布基本保持一致,這是因為在理想模型條件下,單一預報狀態(tài)同樣位于參考狀態(tài)的吸引子上,其大樣本集合等于參考狀態(tài)的吸引子集;而集合平均預報,其概率分布隨時間呈現(xiàn)出值域變窄、峰值變大的特點,這是由于集合平均的非線性濾波作用,使得集合平均預報狀態(tài)的取值縮小到吸引子上的較小區(qū)域內(即吸引子的一個子集上).意味著隨預報時間的推移,單一預報傾向于選擇混沌吸引子上的隨機狀態(tài)進行預報,而集合平均預報則傾向于選擇吸引子子集上的隨機狀態(tài)進行預報.

圖6 Lorenz96模型變量X概率分布,圖中標注為積分步長,黑色虛線為變量X的平均值Fig.6.Probability(%)distribution of variableXof Lorenz96 model with different lengths of integration series.The mean state value(black dashed line)ofXis 2.33.

圖7 Lorenz96模型104個實驗個例對應的變量X的真實狀態(tài)(黑線)、單一預報狀態(tài)(藍線)及NLLV集合平均預報狀態(tài)(紅線)的概率隨時間的變化Fig.7.Probability(%)distribution of the true states(black line,left hand scale),single(blue line,left hand scale)and NLLV ensemble mean(red line,right hand scale)of variableXof Lorenz96 model over the total 104cases as a function of time.

4 結 論

本文利用Lorenz63和Lorenz96理想模型,通過選取吸引子上不同狀態(tài)作為實驗個例,研究了不同實驗個例的集合預報與單一預報的預報表現(xiàn),并分析了集合平均具有較小預報誤差的可能原因,主要研究結論如下.

1)從平均效果看,NLLV和BGM集合平均預報的均方根誤差和型異常相關較單一預報有明顯的改善.預報初期,集合平均與單一預報結果接近;隨著時間的推移,集合平均對預報的改善效果越顯著.

2)不同的實驗個例,其集合預報與單一預報的表現(xiàn)不同.預報早期,集合預報與單一預報的預報技巧相當,隨著時間的推移,集合預報優(yōu)于單一預報的實驗個例數(shù)逐漸增多.在Lorenz63模型中,預報第5 tus時刻,NLLV和BGM集合平均優(yōu)于單一預報的實驗個例數(shù)分別為65%和62%左右;在Lorenz96模型中,預報第12 d時刻,NLLV和BGM集合平均優(yōu)于單一預報的實驗個例數(shù)分別為95%和86%左右.

3)就概率分布(f)而言,單一預報狀態(tài)的f和真實狀態(tài)基本一致,不隨時間變化,而集合平均預報狀態(tài)f隨時間呈現(xiàn)出值域變窄、峰值變大的特點,意味著隨預報時間的延長,單一預報傾向于選擇混沌吸引子上的隨機狀態(tài)進行預報,而集合平均預報傾向于選擇吸引子子集上的隨機狀態(tài)進行預報,這可能是集合平均誤差小于單一預報的原因,但具體的物理機制仍需進一步探討.

上述基于Lorenz模型進行的對比研究結果,對實際的預報業(yè)務啟發(fā)包括:集合預報和單一預報的預報技巧會隨著事件的變化而變化,對不同事件的預報,應該綜合考慮集合預報與單一預報的預報效果;由于集合平均更趨向于系統(tǒng)的平均狀態(tài),對于系統(tǒng)中的極端事件,集合預報的預報效果可能會弱于甚至是差于單一預報.

集合預報作為減小預報結果不確定的可行辦法,值得進一步利用更加復雜的模型(如Weather Research and Forecasting,WRF模式),以不同天氣事件為例,研究其預報表現(xiàn),為集合預報理論發(fā)展和預報技巧的改善提供科學依據(jù).

[1]Lorenz E N 1963J.Atmos.Sci.20 130

[2]Chou J F 1990New Advances of Atmospheric Dynamic(Lanzhou:Lanzhou University Press)p214(in Chinese)[丑紀范1990大氣動力學的新進展(蘭州:蘭州大學出版社)第214頁]

[3]Li J P,Chou J 1996Chin.Sci.Bull.41 587

[4]Epstein E S 1969Tellus21 739

[5]Leith C E 1974Mon.Wea.Rev.102 409

[6]Fritsch J M,Hilliker J,Ross J 2000Wea.Forecasting15 571

[7]Visloscky R L,Fritsch J M 1995Bull.Amer.Meteor.Soc.76 1157

[8]Hoffman R N,Kalnay E 1983Tellus35A 100

[9]Toth Z,Kalnay E 1993Bull.Amer.Meteor.Soc.74 2317

[10]Toth Z,Kalnay E 1997Mon.Wea.Rev.125 3297

[11]Molteni F,Palmer T N 1993Q.J.R.Meteorol.Soc.119 269

[12]Molteni F,Buizza R,Palmer T N,Petroliagis T 1996Q.J.R.Meteorol.Soc.122 73

[13]Buizza R 1996Mon.Wea.Rev.125 99

[14]Houtekamer P L,Lefaivre L,Derome J,Ritchie H,Mitchell H L 1996Mon.Wea.Rev.124 1225

[15]Houtekamer P L,Mitchell H L 1998Mon.Wea.Rev.126 796

[16]Mu M,Jiang Z N 2008Chin.Sci.Bull.53 2062

[17]Duan W S,Mu M 2009Sci.China Ser.D:Earth Sci.52 883

[18]Evensen G 2003Ocean Dyn.53 343

[19]Wang X,and Bishop C 2003J.Atmos.Sci.60 1140

[20]Feng J,Ding R Q,Liu D Q,Li J P 2014J.Atmos.Sci.71 3554

[21]Feng J,Ding R Q,Li J P,Liu D Q 2016Adv.Atmos.Sci.33 1036

[22]Ding R Q,Li J P,Li B S 2017Adv.Atmos.Sci.34 1027

[23]Anderson J L 1997Mon.Wea.Rev.125 2969

[24]Feng G L,Dong W J 2003Acta Phys.Sin.52 2347(in Chinese)[封國林,董文杰 2003物理學報 52 2347]

[25]Bowler N E 2006Tellus58A 538

[26]Li Z C,Chen D H 2002J.Appl.Meteor.Sci.13 1(in Chinese)[李澤椿,陳德輝2002應用氣象學報13 1]

[27]He W P,Feng G L,Dong W J,Li J P 2006Acta Phys.Sin.55 969(in Chinese)[何文平,封國林,董文杰,李建平2006物理學報55 969]

[28]Ma J H,Zhu Y J,Wang P X,Duan M J 2011Trans.Atmos.Sci.34 370(in Chinese)[麻巨慧,朱躍建,王盤興,段明鏗2011大氣科學學報34 370]

[29]Tan N,Chen J,Tian H 2013Meteor.Mon.39 543(in Chinese)[譚寧,陳靜,田華2013氣象 39 543]

[30]Zheng Z H,Feng G L,Huang J P,Chou J F 2012Acta Phys.Sin.61 199203(in Chinese)[鄭志海,封國林,黃建平,丑紀范2012物理學報61 199203]

[31]Zhang H B,Chen J,Zhi X F,Long K J,Wang Y N 2014Trans.Atmos.Sci.37 276(in Chinese)[張涵斌,陳靜,智協(xié)飛,龍柯吉,王亞男2014大氣科學學報37 276]

[32]Ding R Q,Li J P 2009Acta Meteor.Sin.67 241(in Chinese)[丁瑞強,李建平 2009氣象學報 67 241]

[33]Ding R Q,Li J P 2009Acta Meteor.Sin.67 343(in Chinese)[丁瑞強,李建平 2009氣象學報 67 343]

[34]Ding R Q,Li J P,Ha K J 2008J.Geophys.Res.113 D24112

[35]Li J P,Ding R Q 2011Mon.Wea.Rev.139 3265

[36]Lorenz E N,Emanuel K A 1998J.Atmos.Sci.55 399

[37]Lorenz E N 1995Seminar on PredictabilityShinf i eld Park,Reading,United Kingdom,September 4–8,1995 p1

[38]Buizza R,Houtekamer P L,Toth Z,Pellerin G,Wei M,Zhu Y 2005Mon.Wea.Rev.133 1076

[39]Houtekamer P L,Derome J 1994Mon.Wea.Rev.122 2179

[40]Zhuang Z R,Xue J S,Li X L 2011Acta Meteor.Sin.69 620(in Chinese)[莊照榮,薛紀善,李興良2011氣象學報69 620]

[41]Eckmann J P,Ruelle D 1985Rev.Mod.Phys.57 617