康俊鋒 馮松江 鄒倩 李艷杰 丁瑞強 鐘權(quán)加?
1) (江西理工大學(xué)土木與測繪工程學(xué)院,贛州 341000)
2) (中國科學(xué)院大氣物理研究所,大氣科學(xué)和地球流體力學(xué)數(shù)值模擬國家重點實驗室,北京 100029)
3) (北京師范大學(xué),地表過程與資源生態(tài)國家重點實驗室,北京 100875)
4) (中國科學(xué)院大學(xué)地球科學(xué)學(xué)院,北京 100049)
基于Lorenz96 模型初步探討了機器學(xué)習(xí)算法提高非線性局部Lyapunov 向量(NLLV)集合預(yù)報效果的可行性和有效性.結(jié)果表明:基于嶺回歸算法和NLLV 集合預(yù)報結(jié)果建立的機器學(xué)習(xí)模型(Ens-ML)能夠有效提高整體預(yù)報技巧,而且優(yōu)于集合平均預(yù)報(EnsAve)、控制預(yù)報(Ctrl)以及基于Ctrl 結(jié)果建立的機器學(xué)習(xí)模型(Ctrl-ML).同時,還發(fā)現(xiàn)Ens-ML 的預(yù)報技巧改進(jìn)程度依賴于集合成員的數(shù)量,即增加集合成員數(shù)有助于提高Ens-ML 模型的整體預(yù)報準(zhǔn)確率.通過對比個例預(yù)報表現(xiàn)得到,隨著預(yù)報時間延長,Ens-ML,Ctrl-ML和EnsAve 的個例預(yù)報誤差逐漸小于Ctrl.進(jìn)一步分析Ens-ML,Ctrl-ML 和EnsAve 預(yù)報的吸引子,發(fā)現(xiàn)它們的概率分布的值域收縮、峰度增大并向平均值靠攏,尤其Ens-ML 的表現(xiàn)更為明顯.
近幾十年來,得益于計算機硬件的快速發(fā)展、數(shù)值模式和資料同化技術(shù)的研發(fā)和應(yīng)用,數(shù)值預(yù)報取得了長足的進(jìn)步,提升了天氣預(yù)報和氣候預(yù)測的準(zhǔn)確率.然而,大氣作為一個具有混沌特性的非線性系統(tǒng),其預(yù)報結(jié)果存在一定的不確定性[1].最早,Leith[2]和Epstein[3]提出了集合預(yù)報的思想和方法來應(yīng)對數(shù)值預(yù)報中的不確定性問題,其基本思路是在初始狀態(tài)上疊加一組擾動生成初始集合,并分別積分各個初始成員得到一組預(yù)報結(jié)果的集合.一方面,集合平均的非線性濾波作用可以過濾不可預(yù)報的噪音,保留可預(yù)報的有用信息;另一方面,根據(jù)集合成員的預(yù)報結(jié)果可獲得概率預(yù)報.自20 世紀(jì)90 年代以來,集合預(yù)報在天氣預(yù)報與氣候預(yù)測領(lǐng)域得到了廣泛的研究和應(yīng)用,眾多研究和業(yè)務(wù)應(yīng)用的結(jié)果都表明集合預(yù)報是減小預(yù)報不確定性、提高預(yù)報水平的有效途徑[4?6].
針對數(shù)值預(yù)報中的初值不確定性和模式本身的不確定性,前人相應(yīng)發(fā)展了初值擾動和模式物理擾動以及多模式集合預(yù)報方法[7].假設(shè)在模式完美情況下,初始集合生成的關(guān)鍵是如何產(chǎn)生能夠很好描述初始狀態(tài)不確定性的擾動場.為此,國內(nèi)外學(xué)者相繼發(fā)展了多種初值集合生成方法,代表性的包括繁殖向量(BV)[8,9]、奇異向量(SV)[10]和集合卡爾曼濾波(EnKF)[11,12]方法.其中,BV 和SV 方法發(fā)展較為成熟,它們都是基于誤差增長動力學(xué)理論發(fā)展起來的方法,主要通過抓住增長的擾動結(jié)構(gòu)和演變來反映數(shù)值預(yù)報誤差的不確定性.盡管BV和SV 方法在業(yè)務(wù)上得到較好的應(yīng)用,但仍存在一定的不足.例如,SV 方法基于線性誤差增長理論,具有線性約束的局限性,不能很好地描述誤差非線性增長特征.為此,Duan 和Huo[13]發(fā)展了CNOP 方法來表征誤差的非線性增長模態(tài).與SV 和CNOP方法相比,BV 方法無需復(fù)雜的切線性和伴隨模式,計算非常簡便和省時.但是BV 擾動是在動力系統(tǒng)中相同的基流演化得到,導(dǎo)致產(chǎn)生的初始擾動結(jié)構(gòu)較為相似,獨立性不足[14].考慮BV 方法的優(yōu)點和不足,最近我國學(xué)者提出了非線性局部Lyapunov向量(NLLV)[15,16]集合預(yù)報新方法,該方法利用和BV 方法類似的簡便的繁殖思想生成擾動,能夠很好地抓住非線性動力系統(tǒng)中的增長誤差結(jié)構(gòu),同時又通過正交化彌補BV 方法產(chǎn)生的擾動結(jié)構(gòu)依賴性較強的不足,有助于進(jìn)一步提升集合初值擾動的產(chǎn)生質(zhì)量.
近年來,隨著大數(shù)據(jù)科學(xué)的發(fā)展,越來越多的研究將機器學(xué)習(xí)算法和人工智能技術(shù)應(yīng)用于挖掘海量的氣象數(shù)據(jù)信息并加以利用,以提高天氣預(yù)報和氣候預(yù)測的準(zhǔn)確率[17?19].一方面,既有利用單一數(shù)值模式的不同氣象要素構(gòu)建預(yù)報預(yù)測模型,例如,利用機器學(xué)習(xí)中的隨機森林算法構(gòu)建短時地面風(fēng)場預(yù)報模型,能夠有效預(yù)報地面1—6 h 風(fēng)場變化[20];另一方面,也有基于單個數(shù)值模式集合預(yù)報結(jié)果或者多模式的集合預(yù)報結(jié)果,結(jié)合機器學(xué)習(xí)算法建立集合預(yù)報機器學(xué)習(xí)模型,從而達(dá)到提高預(yù)報預(yù)測準(zhǔn)確率的目的.比如,通過機器學(xué)習(xí)算法對華北氣溫的多模式集合預(yù)報結(jié)果進(jìn)行建模訂正,其訂正效果明顯優(yōu)于數(shù)值模式的單一預(yù)報和多模式的集合平均[21].此外,也有學(xué)者基于深度學(xué)習(xí)模型構(gòu)建了長期氣候預(yù)測模型,實現(xiàn)了對ENSO 的有效預(yù)測[22].
鑒于NLLV 是一種新的集合預(yù)報方法,具有一定的理論創(chuàng)新性,同時機器學(xué)習(xí)方法又在氣象領(lǐng)域具有廣闊的應(yīng)用前景.因此,本文將基于簡單的Lorenz96 模型開展集合預(yù)報試驗,然后結(jié)合機器學(xué)習(xí)算法分別利用NLLV 的集合預(yù)報和單一控制預(yù)報建立機器學(xué)習(xí)模型,并從整體預(yù)報效果、不同個例預(yù)報誤差大小、吸引子概率分布以及集合成員個數(shù)的使用等多個不同角度探討機器學(xué)習(xí)算法和NLLV 集合預(yù)報相結(jié)合的可行性和有效性.
Lorenz96 模型具有對初值極端敏感性和非線性的特點,能夠作為大氣系統(tǒng)的低階近似,廣泛應(yīng)用于大氣可預(yù)報性和集合預(yù)報研究[15,23],該模型動力方程組可表示為

其中,Xi(i=1,2,···,40)為模型狀態(tài)變量,定常強迫項F=8,數(shù)值積分采用四階龍格-庫塔格式,積分步長為0.05 時間單位(tus).在Lorenz96 模型的混沌吸引子上每間隔0.05 tus 選取一個狀態(tài)作為集合預(yù)報試驗的初始態(tài),本研究共選取了10000個初始態(tài)作為集合預(yù)報個例,每個個例的積分步數(shù)為41.在Lorenz96 模型中每間隔0.2 tus 大致對應(yīng)1 d[24].因此,積分41 步共為2.05 tus,大致對應(yīng)10 d.
NLLV 是一種生成集合預(yù)報初始擾動的新方法[15,16],該方法產(chǎn)生的初始集合成員NLLVs 具有正交性、獨立性以及流依賴性等特點.因此,它們既能夠描述不同方向的誤差增長率,還能夠很好地反映誤差結(jié)構(gòu)隨天氣、氣候等混沌系統(tǒng)時空演化的特征[25].基于上述優(yōu)勢,NLLV 方法逐漸應(yīng)用于Lorenz 模型、準(zhǔn)地轉(zhuǎn)正壓模式、中等復(fù)雜程度的ENSO 動力耦合模式(Zebiak-Cane 模式)以及新一代中尺度數(shù)值模式(WRF 模式)的集合預(yù)報研究,并發(fā)現(xiàn)其能夠有效地提高預(yù)報技巧[15,26,27].NLLV 方法產(chǎn)生擾動制作集合預(yù)報的流程如圖1所示,當(dāng)NLLV 繁殖開始時,將一組初始隨機擾動添加到分析狀態(tài)上進(jìn)行積分,再對積分的結(jié)果進(jìn)行正交化處理并尺度化至初始擾動大小,然后又將其疊加到新的分析狀態(tài)上,經(jīng)過多個繁殖循環(huán)得到NLLVs.參考Feng 等[15]的研究,當(dāng)集合成員達(dá)到一定數(shù)量時,使用更多的成員也只能有限地提高預(yù)報技巧,因此本研究使用NLLV 方法共生成6 個集合成員.

圖1 NLLVs 擾動生成示意圖[15]Fig.1.Schematic diagram of the generation of NLLVs[15] .
嶺回歸屬于多元線性回歸方法的一種,其在最小二乘法的基礎(chǔ)上改進(jìn)而來,能夠很好地解決多重共線性數(shù)據(jù)集的回歸問題[28,29].考慮到本研究的NLLV 集合預(yù)報試驗數(shù)據(jù)集中各預(yù)報成員之間呈現(xiàn)強相關(guān)性,具有顯著的共線性特征(圖2),所以將該方法作為構(gòu)建機器學(xué)習(xí)模型的算法.在嶺回歸模型的構(gòu)建過程中將預(yù)報成員作為特征變量xi(i=1,2,3,···,7),其模型框架可表示為

圖2 各集合成員間的相關(guān)系數(shù)矩陣Fig.2.Correlation coefficient matrix of ensemble members.

即

模型的損失函數(shù)為

其中wi為特征變量的回歸系數(shù),b為偏置,y為真值,λ為正則項系數(shù);w和x分別是wi和xi的矩陣.用梯度下降算法對損失函數(shù)求最小解可得w和b,確定嶺回歸模型.
圖3 給出了機器學(xué)習(xí)的建模流程.首先,對集合成員進(jìn)行標(biāo)準(zhǔn)化預(yù)處理;然后,將集合預(yù)報試驗的控制預(yù)報成員(Control,Ctrl)作為特征變量,真值作為目標(biāo)變量,建立基于控制預(yù)報的機器學(xué)習(xí)模型(Control machine learning,Ctrl-ML);接著,將集合預(yù)報試驗的集合成員作為特征變量,真值作為目標(biāo)變量,建立基于集合成員的機器學(xué)習(xí)模型(Ensemble machine learning,Ens-ML);上述建模過程還需要分別將兩個模型的數(shù)據(jù)集劃分為訓(xùn)練集和測試集,其中訓(xùn)練集占70% (7000 個),測試集占30% (3000 個).接著,分別對上述兩個模型進(jìn)行訓(xùn)練、測試,在訓(xùn)練模型時將其迭代次數(shù)參數(shù)“max_iter”設(shè)置為None,模型將自動尋得最優(yōu)迭代次數(shù),正則項系數(shù)“alpha”為1.0.最后,對模型的測試結(jié)果進(jìn)行評估并輸出,評估方法見2.4 節(jié).

圖3 機器學(xué)習(xí)模型構(gòu)建流程圖Fig.3.Process of machine learning.
通過對所有集合成員做算術(shù)平均,得到集合平均預(yù)報(ensemble average,EnsAve).即

其中,Vi代表集合成員的值,N代表集合成員個數(shù).
此外,為了從不同角度來綜合評估不同方法的預(yù)報表現(xiàn),本文主要使用了可決系數(shù)(R2)、平均絕對誤差(MAE)、均方誤差(MSE)、均方根誤差(RMSE)和型異常相關(guān)系數(shù)(PAC)等預(yù)報評估方法[30,31].
圖4 給出了一組不同預(yù)報方法下Lorenz96 模型X變量的預(yù)報時間序列,可以發(fā)現(xiàn)4 種方法在不同的預(yù)報時段各有優(yōu)劣.在預(yù)報初始時刻0—5 tus 時,Ctrl 和EnsAve 的效果優(yōu)于Ctrl-ML 和Ens-ML;在預(yù)報中期20—28 tus,Ctrl 和Ctrl-ML 的預(yù)報值最接近真值,Ens-ML 和EnsAve 則略遜;但是在31 tus 時的波峰處,Ctrl 試驗的結(jié)果卻是誤差最大的;到后期35—40 tus 時,EnsAve和Ens-ML 的預(yù)報結(jié)果明顯優(yōu)于Ctrl 和Ctrl-ML.通過R2和MAE 等評價指標(biāo)從整體表現(xiàn)來看,Ens-ML的R2最大,MAE 和RMSE 最小,表明Ens-ML 的整體預(yù)報技巧是最高的,EnsAve 次之,Ctrl 的預(yù)報技巧最差(見表1).

圖4 不同預(yù)報方法下Lorenz96 模型的X 變量時間序列 (a) Ctrl;(b) EnsAve;(c) Ctrl-ML;(d) Ens-ML(黑線為真值的時間序列)Fig.4.The time series of X variable of Lorenz96 model for different forecast methods:(a) Ctrl;(b) EnsAve;(c) Ctrl-ML;(d) Ens-ML.Black line represents time series of true values.

表1 不同預(yù)報方法的預(yù)報結(jié)果比較Table 1.Evaluation of forecast results in different forecasting methods.
如圖5(a)所示,由于兩個機器學(xué)習(xí)模型均在預(yù)報初期出現(xiàn)過擬合現(xiàn)象,從而導(dǎo)致其平均RMSE大于Ctrl 和EnsAve 預(yù)報.然而,隨著預(yù)報時間的延長,Ctrl-ML 和Ens-ML 的RMSE 逐漸分別小于Ctrl 和EnsAve 的誤差,尤其是Ens-ML 表現(xiàn)最好.此外,根據(jù)預(yù)報序列與真值的PAC 可以發(fā)現(xiàn),Ens-ML 優(yōu) 于EnsAve,而Ctrl-ML 與Ctrl 的PAC 曲線重合,這可能是在Ctrl-ML 機器學(xué)習(xí)模型中,只有Ctrl 作為其特征變量,Ctrl-ML 只能學(xué)習(xí)到Ctrl 的模式信息,導(dǎo)致它們的PAC 與真值的相關(guān)性較為一致.但是,Ens-ML 則能綜合多個集合成員的信息,使得最后效果優(yōu)于EnsAve (圖5(b)).進(jìn)一步比較3 種不同方法相對于Ctrl 的改進(jìn)程度可以發(fā)現(xiàn),雖然兩個機器學(xué)習(xí)模型在預(yù)報初期預(yù)報誤差偏大,提升比例為負(fù)(圖6),但是,隨著時間的演變EnsAve,Ctrl-ML 和Ens-ML 的預(yù)報誤差逐漸小于Ctrl,改進(jìn)程度逐漸提高,直到最后Ens-ML 提升比例最高可達(dá)3%以上,EnsAve 次之,Ctrl-ML 提升最小.上述結(jié)果表明,機器學(xué)習(xí)對集合預(yù)報和控制預(yù)報的整體性能是有提高作用的,并且以集合成員為特征建立的機器學(xué)習(xí)模型明顯優(yōu)于基于單一控制預(yù)報試驗建立的機器學(xué)習(xí)模型.

圖5 不同預(yù)報方法得到的平均RMSE (a)和平均PAC (b)隨時間步長的變化Fig.5.The average RMSE (a) and average PAC (b) for different forecast methods.

圖6 EnsAve,Ens-ML,Ctrl-ML 相對于Ctrl 預(yù)報的改進(jìn)程度Fig.6.The improvement of the EnsAve,Ens-ML and Ctrl-ML compared to Ctrl.
圖7 給出了不同預(yù)報方法的誤差概率分布,可以看出Ens-ML 相比EnsAve 略微提高了誤差在[0,2.5]區(qū)間的分布概率,降低了在[2.5,4.5]區(qū)間的概率,在高誤差區(qū)間[4.5,6]與EnsAve 基本保持一致.Ctrl-ML 相比Ctrl 則更加明顯,顯著提高了[1,3]區(qū)間的誤差概率,降低了[3,6]區(qū)間的誤差概率.此外,根據(jù)圖8 可以發(fā)現(xiàn)機器學(xué)習(xí)模型的預(yù)報效果與集合成員的數(shù)量有重要聯(lián)系.在Ens-ML 模型中,預(yù)報結(jié)果與真值的R2隨著集合成員的增加而增大.這意味著機器學(xué)習(xí)依賴NLLV 的集合成員,集合成員的增加能夠顯著提高模型的預(yù)報準(zhǔn)確度.這一結(jié)果也間接解釋了Ctrl-ML 模型預(yù)報效果明顯不如Ens-ML 的原因.相對于Ctrl試驗單一成員,集合成員更能夠反映出大氣系統(tǒng)的不確定性狀態(tài).集合成員的增加能夠給機器學(xué)習(xí)模型提供更多的大氣特征信息,從而提高模型的預(yù)報精度.綜上所述,整體來看機器學(xué)習(xí)模型提高了誤差在低值區(qū)的概率,降低了高值區(qū)的概率,有助于提升整體預(yù)報效果,并且集合成員數(shù)量對機器學(xué)習(xí)模型的效果起著至關(guān)重要的影響.

圖7 不同預(yù)報方法的誤差概率分布Fig.7.The probability distribution of forecast errors for different methods.

圖8 Ens-ML 模型的R2 隨集合成員數(shù)的變化Fig.8.Changes of the R2 with the number of ensemble member used in the Ens-ML model.
前人的研究曾指出集合平均預(yù)報整體預(yù)報技巧優(yōu)于單一控制預(yù)報,但也存在一定比例的個例其單一控制預(yù)報優(yōu)于集合預(yù)報[23].下文將進(jìn)一步比較EnsAve,Ens-ML,Ctrl-ML 模型與Ctrl 單一控制預(yù)報在3000 個測試個例中的預(yù)報誤差大小.
在試驗初期0 tus 時,EnsAve 和Ens-ML 的試驗個例與Ctrl 的試驗個例主要集中在對角線附近,說明初期時這3 種方法的個例誤差大小相當(dāng)(圖9(a),(b),(c)).隨著預(yù)報時間的延長,大多數(shù)個例的EnsAve 和Ens-ML 模型誤差小于Ctrl,該結(jié)果與前人研究一致.除此之外,由于Ctrl-ML 模型的構(gòu)建過程使用的特征變量只有Ctrl 成員,使得機器學(xué)習(xí)算法在對其訓(xùn)練過程中只有單一特征信息與真值對應(yīng)學(xué)習(xí),導(dǎo)致與Ctrl 自身個例的對比中整體呈現(xiàn)出線性變化趨勢.

圖9 試驗個例在不同時刻的EnsAve,Ctrl-ML 和Ens-ML 與Ctrl 的預(yù)報誤差 (a)—(c) 0 tus;(d)—(f) 10 tus;(g)—(i) 20 tus;(j)—(l) 30 tus;(m)—(o) 40 tusFig.9.Scatterplot of forecast error at different leading times between the EnsAve,Ctrl-ML,Ens-ML and the Ctrl,respectively:(a)–(c) 0 tus;(d)–(f) 10 tus;(g)–(i) 20 tus;(j)–(l) 30 tus;(m)–(o) 40 tus .
Lorenz96 模型中吸引子的概率分布可以揭示系統(tǒng)變量在演化過程中的狀態(tài)特征,圖10 給出了3000 測試個例采用4 種不同預(yù)報方法在整個預(yù)報周期中Lorenz96 模型的X變量的概率分布.由此可以看出,Ens-ML,Ctrl-ML 的概率分布趨向與真值是一致的,這表明本試驗建立的機器學(xué)習(xí)模型是可靠的,吸引子的系統(tǒng)演化軌跡是基本不變的.隨著預(yù)報時間的延長,不同時刻Lorenz96 模型X變量的概率分布總體上分布一致,但EnsAve,Ens-ML 和Ctrl-ML 三種方法的概率分布逐漸向中間均值部分收縮,左右分布概率降低,值域變窄,峰度值越來越大,Ctrl 則始終與真值基本保持一致.其中Ens-ML 和EnsAve 的值域收縮最為明顯,Ctrl-ML 相較于Ctrl 收縮程度大.這說明機器學(xué)習(xí)模型在訓(xùn)練學(xué)習(xí)過程中,有著同集合預(yù)報一樣偏向平均態(tài)的特征且比EnsAve 效果更好,從而提高了預(yù)報技巧,但左右極值分布概率的降低意味對極端事件的預(yù)報技巧是下降的.

圖10 不同預(yù)報方法中X 變量狀態(tài)的概率分布隨時間的變化Fig.10.Probability distributions of X variables for different leading times.
圖11 給出了4 種方法的預(yù)報值與真值的分布對比,也可以看出Ens-ML 和EnsAve 更集中更靠近真值分布,Ctrl-ML 次之,而Ctrl 較前三者分布最為發(fā)散,預(yù)報技巧最差.可以發(fā)現(xiàn)在極值兩端,更多Ctrl 預(yù)報點偏向極值,而EnsAve 和Ens-ML則明顯比Ctrl 要小得多,這說明EnsAve 和Ens-ML與真值的誤差更小,比Ctrl 更接近真值.

圖11 不同預(yù)報方法的結(jié)果與真實狀態(tài)的對比分布Fig.11.Scatterplot of the forecast value.
綜上所述,隨著預(yù)報時間的推移,多數(shù)個例的Ctrl 誤差逐漸大于EnsAve 和機器學(xué)習(xí)模型,并且機器學(xué)習(xí)模型和EnsAve 的吸引子的概率分布出現(xiàn)值域變窄,峰度值增大的特點,總得說來Ens-ML,Ctrl-ML 和EnsAve 的預(yù)報技巧均優(yōu)于Ctrl,其中又以Ens-ML 預(yù)報效果最佳.
本文基于Lorenz96 模型,借助嶺回歸算法構(gòu)建了基于控制預(yù)報和NLLV 集合預(yù)報的機器學(xué)習(xí)模型,初步探究了機器學(xué)習(xí)算法改進(jìn)NLLV 集合預(yù)報效果的可行性,得到的主要結(jié)論如下.
1) 機器學(xué)習(xí)能夠有效提高NLLV 集合預(yù)報的預(yù)報技巧,且其預(yù)報改進(jìn)程度明顯依賴于構(gòu)建機器學(xué)習(xí)模型的集合成員數(shù),增加集合成員數(shù)能夠減小預(yù)報初期的過擬合現(xiàn)象,提高預(yù)報結(jié)果與真值之間的相關(guān)系數(shù),究其原因可能是較多的集合成員能夠提供更多的特征變量信息,使得機器學(xué)習(xí)訓(xùn)練的模型精度更高.
2) 就試驗個例表現(xiàn)而言,隨著時間的演變,控制預(yù)報誤差大于其他3 種方法的個例數(shù)增多,在試驗中后期尤為明顯.從吸引子的概率分布來看,機器學(xué)習(xí)和集合平均都具有值域收縮、峰度增大并向平均值靠攏的特點,且機器學(xué)習(xí)比集合平均表現(xiàn)更突出,說明機器學(xué)習(xí)模型的預(yù)報結(jié)果優(yōu)于集合平均和控制預(yù)報.
本研究基于Lorenz96 模型利用機器學(xué)習(xí)算法訂正NLLV 集合預(yù)報結(jié)果,期望進(jìn)一步改進(jìn)其預(yù)報技巧.研究發(fā)現(xiàn),機器學(xué)習(xí)模型在中后期對集合預(yù)報有較為明顯的改進(jìn),同時在預(yù)報初期有一定的過擬合現(xiàn)象.雖然,增加集合成員數(shù)目能夠一定程度上解決初期的過擬合現(xiàn)象,但也意味著會增加計算量,需要更多的計算資源.