李孝涌, 陳科藝, 李熙晨
(1.成都信息工程大學(xué)大氣科學(xué)學(xué)院,四川 成都 610225;2.中國科學(xué)院大氣物理研究所,北京 100029)
厄爾尼諾與南方濤動(El Nino-southern oscillation,ENSO)是發(fā)生在熱帶中東太平洋的海溫、風(fēng)場異常現(xiàn)象[1],是熱帶太平洋主要的氣候變率,對全球的氣候有深遠的影響。海表溫度異常作為延伸期重要的下墊面信號,能通過大氣遙相關(guān)作用對中國的氣溫和降水造成影響[2],因此,對ENSO的準(zhǔn)確預(yù)測能提高中長期氣候的預(yù)測技巧。
目前,ENSO的預(yù)報方法主要分為動力方法和統(tǒng)計方法兩大類。動力方法采用動力學(xué)方程,對海洋、大氣、陸地等圈層及其相互作用進行建模,并利用計算機逐步積分,模擬大氣的演變,計算出未來的氣候狀況。歐洲中期天氣預(yù)報中心(ECMWF)于2017年推出的第五代季節(jié)預(yù)報系統(tǒng)SEAS5包含了大氣、海洋和海冰等部分,是一個耦合動力模型。對SEAS5的評估顯示,其在ENSO預(yù)報方面繼續(xù)顯示出特別的優(yōu)勢。在預(yù)報時效為9個月時,SEAS5對ENSO預(yù)報的相關(guān)系數(shù)在0.7以上[3]。相對而言,統(tǒng)計學(xué)模式通過一系列統(tǒng)計學(xué)方法,如線性轉(zhuǎn)置模型(LIM)、非線性典型相關(guān)分析(NLCCA)、奇異譜分析(SSA)等[4]來預(yù)報未來的氣候狀況。Xue等[5]利用線性Markov模型,采用海表溫度、海平面高度和風(fēng)應(yīng)力作為預(yù)報因子,建立了預(yù)報模型。在預(yù)報時效為6個月時,其預(yù)報相關(guān)技巧為0.8。Kondrashov[6]采用多層次回歸分析來獲得ENSO的隨機強迫模型。在預(yù)報時效為6個月時,其相關(guān)系數(shù)超過0.6。
近些年,隨著機器學(xué)習(xí)技術(shù)的快速發(fā)展和計算機硬件計算速度的不斷提高,基于神經(jīng)網(wǎng)絡(luò)的預(yù)測方法被廣泛應(yīng)用在各個領(lǐng)域,神經(jīng)網(wǎng)絡(luò)在氣象領(lǐng)域中的應(yīng)用也正日益受到關(guān)注。ENSO作為海洋-大氣相互作用的產(chǎn)物,具有明顯的非線性特征。而神經(jīng)網(wǎng)絡(luò)中引入的非線性的激活函數(shù),使其具備了擬合任意函數(shù)的能力[7],從而能夠完成非線性回歸的任務(wù)。
前人在將機器學(xué)習(xí)用于ENSO預(yù)報上已做了諸多有益的嘗試。許柏寧等[8]使用序列到序列(Seq2Seq)模型來預(yù)測區(qū)域SSTA,結(jié)果在提前7個月以上的預(yù)報上,均方根誤差(RMSE)比動力學(xué)模式小,但在峰值處誤差較大。Dandan He等[9]同樣使用了基于卷積長短時記憶神經(jīng)網(wǎng)絡(luò)(ConvLSTM)的Seq2Seq模型用于預(yù)報Nino3.4指數(shù),其效果比長短期記憶網(wǎng)絡(luò)(LSTM)、動力學(xué)模式的決定性預(yù)報和集合預(yù)報好。Clifford Broni Bedaiko等[10]使用“氣候復(fù)雜網(wǎng)絡(luò)”(Climate Complex Network)提取氣象數(shù)據(jù)特征,用提取得到的特征作為預(yù)報因子,使用LSTM預(yù)報Nino3.4指數(shù)。Ham等[11]以過去3個月的SST和海洋熱含量(heat content,HC)為預(yù)測因子,訓(xùn)練了一個CNN,該CNN在Nino3.4指數(shù)與厄爾尼諾類型的預(yù)測上明顯優(yōu)于動力學(xué)模型,可提前18個月對ENSO做出預(yù)測。
目前,機器學(xué)習(xí)用于ENSO預(yù)報時,與集合預(yù)報方法結(jié)合比較少。同時各種因素對機器學(xué)習(xí)預(yù)報效果的影響也有待探討?;谶@個原因,文中建立了以深度卷積神經(jīng)網(wǎng)絡(luò)(CNN)為基礎(chǔ)的ENSO預(yù)報模型,通過對照實驗確定最佳的模型,并討論模型在后報實驗中給出的結(jié)果。
用于訓(xùn)練和測試神經(jīng)網(wǎng)絡(luò)的數(shù)據(jù)包括觀測資料和模式資料。其中,觀測資料來自美國國家環(huán)境預(yù)報中心(NCEP)提供的全球海洋資料同化系統(tǒng)(global ocean data assimilation system,GODAS)產(chǎn)品。GODAS產(chǎn)品包含海水鹽度、海水溫度、洋流和海平面高度等多種變量,范圍涵蓋75°S ~65°N,最大分辨率為1°×1°[12]。模式資料來自第六次國際耦合模式比較計劃(CMIP6)。Mckenna等[13]從季節(jié)性、強度和頻率考察了CMIP6模式對ENSO的模擬能力,結(jié)果顯示部分模式對ENSO的模擬在季節(jié)鎖相特征和衰亡階段存在偏差,多數(shù)模式模擬的ENSO頻率與觀測到的相同,說明CMIP6對ENSO有一定的模擬能力。選用的CMIP6模式包括:AWI-CM-1-1-MR、BCC-CSM2-MR、BCC-ESM1、CanESM5、CESM2、CESM2-WACCM、MCM-UA-1-0、MIROC6、NESM3、SAM0-UNICON。不同CMIP6模式之間的水平和垂直分辨率存在差異,因此使用前需統(tǒng)一進行資料預(yù)處理。對于上述所有資料,選取的時間分辨率為逐月數(shù)據(jù),選取的資料變量為海表溫度以及海水位溫。
用于比較預(yù)報效果的ENSO歷史預(yù)報數(shù)據(jù)由哥倫比亞大學(xué)的國際氣候與社會研究所提供,包括了各個動力模式和統(tǒng)計模式的Nino 3.4區(qū)歷史預(yù)報數(shù)據(jù),時間跨度為2002-2019年。模式包含了動力模式和統(tǒng)計模式。模式的具體情況[22]可參照表1及表2。在診斷海-氣相互作用強度時,Nino 3.4指數(shù)由NOAA物理科學(xué)實驗室(PSL)提供,南方濤動指數(shù)(southern oscillation index,SOI)由東英吉利亞大學(xué)氣候研究組提供。
表1 用于比較預(yù)報效果的動力模式
表2 用于比較預(yù)報效果的統(tǒng)計模式
對用于訓(xùn)練神經(jīng)網(wǎng)絡(luò)的模式和觀測資料,選取的空間范圍為中低緯度(60°S~60°N)區(qū)域。首先,將資料統(tǒng)一雙線性內(nèi)插到分辨率為2.5°×2.5°(格點數(shù)114×49)的格點上,之后計算資料的距平值,然后計算3個月的滑動平均值,最后得到的資料為3個月滑動平均距平值。
海表溫度(SST)定義為海洋表層的海水溫度。熱含量(HC)定義為0~300 m的海水位溫積分[14]。Nino 3.4指數(shù)是ENSO的指標(biāo),定義為Nino 3.4區(qū)(170°W ~120°W,5°S~5°N)的緯度加權(quán)平均海溫的距平值[15]。
使用基于卷積神經(jīng)網(wǎng)絡(luò)(CNN)的深度學(xué)習(xí)架構(gòu)。CNN作為神經(jīng)網(wǎng)絡(luò)的一種類型,能有效提取圖像中包含的特征,因此被廣泛應(yīng)用在涉及圖像處理(如圖像識別、目標(biāo)檢測等)的領(lǐng)域[16]。對氣象數(shù)據(jù)來說,可以把某時刻某要素的分布場視作一張圖像,將其作為CNN的輸入。用CNN來解決,實際上就是全球海洋要素場與未來幾個月Nino 3.4區(qū)海溫的非線性回歸問題。
CNN中包含卷積運算和下采樣運算。卷積運算可表示為
另定義代價函數(shù):
其中,C為代價函數(shù)的損失值,是實際值與預(yù)測值之差,是權(quán)重w和偏置b的函數(shù)。優(yōu)化目標(biāo)是尋找代價函數(shù)的最小值,即用梯度下降法尋找使C最小的w,b。
其中η為學(xué)習(xí)率。訓(xùn)練神經(jīng)網(wǎng)絡(luò)時,利用式(3)、(4)可迭代更新權(quán)重w和偏置b。
所用的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)包含2個輸入和9個輸出。輸入1為近3個月全球中低緯度區(qū)域(60°S~60°N)的海表溫度(SST)距平和海洋0~300m熱含量(HC)距平;輸入2為近3個月的Nino 3.4指數(shù);網(wǎng)絡(luò)的輸出為未來9個月Nino 3.4指數(shù)的預(yù)測值。
ENSO是海氣相互作用的結(jié)果,但相比于大氣,海洋存在巨大的熱慣性,運動變化有緩慢性和持續(xù)性的特征,同時其本身也能反映一定大氣環(huán)流變化的信息。Ham等[11]的實驗證明僅使用海洋的數(shù)據(jù)是可行的?;诖?文中同樣使用海洋的數(shù)據(jù)作為輸入。預(yù)測Nino 3.4指數(shù)時,有必要特別考慮近3個月的Nino 3.4指數(shù)的變化。盡管輸入1(中低緯SSTA)包含了Nino 3.4海溫的信息,仍將Nino3.4海溫單獨作為另一個輸入,這樣可以使模型更好地利用過往Nino 3.4指數(shù)的信息。另外,考慮到多數(shù)統(tǒng)計模式都采用近3個月的氣象要素作為輸入,為了方便比較在同一情況下的效果,文中神經(jīng)網(wǎng)絡(luò)的輸入同樣為近3個月的資料。
神經(jīng)網(wǎng)絡(luò)的詳細結(jié)構(gòu)參見表3。其中,v為變量個數(shù),在對照實驗中,有只輸入SST和同時輸入SST和HC兩種情況,因此v的取值為1或2。t為輸入時間步數(shù),實驗中取3個月。
表3 采用的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)
神經(jīng)網(wǎng)絡(luò)中的卷積層起到提取特征的作用,下采樣處理則能通過取最大值保留卷積層提取出來的最明顯的特征。通過卷積運算和下采樣運算,得到的特征圖的尺寸被不斷縮小。相對而言,卷積核的感知范圍不斷擴大。因此,采用卷積神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)能夠從全球氣象要素場中提取小尺度到大尺度的各個不同尺度的特征。
為研究不同算法、資料和變量對模型預(yù)測效果的影響,設(shè)計了如下的對照試驗,除對照變量外,其余條件保持一致。在下述條件下訓(xùn)練神經(jīng)網(wǎng)絡(luò),當(dāng)代價函數(shù)給出的損失值C不再有明顯下降時,則認(rèn)為神經(jīng)網(wǎng)絡(luò)模型收斂,訓(xùn)練已完成。此時終止訓(xùn)練,并將模型用于后報實驗中以檢測模型的效果。
對照實驗一分為兩組,區(qū)別在于使用算法的不同。首先在保持訓(xùn)練集一致的情況下,引入50個模型分別訓(xùn)練。由于神經(jīng)網(wǎng)絡(luò)的初始權(quán)重是隨機賦值的,最終得到的50個模型略有不同。該實驗的其中一組在做預(yù)測時,依次使用每一個網(wǎng)絡(luò)單獨做預(yù)測(下稱單一模型),共預(yù)測50次。評價單一模型時,最終結(jié)果取這50次預(yù)測效果的平均。該實驗的另一組則共同將這50個模型輸出的平均作為結(jié)果,使用多個模型做預(yù)測的做法在機器學(xué)習(xí)領(lǐng)域又被稱為集成學(xué)習(xí),在時間序列預(yù)測上應(yīng)用廣泛,前人的研究顯示采用這種方法對改進預(yù)報結(jié)果能起到一定效果[17-18]。集成學(xué)習(xí)與數(shù)值預(yù)報中的集合預(yù)報的思路相似,不同之處在于集合預(yù)報中的各子成員的初始擾動需要用各種方法專門確定,而集成學(xué)習(xí)是通過訓(xùn)練時引入隨機性來產(chǎn)生不同的子成員。下文將采用集成學(xué)習(xí)的模型稱為集合模型。
對照實驗二分為三組,區(qū)別在于使用資料不同。其中組一采用實際觀測資料,組二采用模式資料,而組三同時采用觀測資料和模式資料。當(dāng)同時使用觀測資料和模式資料訓(xùn)練網(wǎng)絡(luò)時,采取了機器學(xué)習(xí)領(lǐng)域中遷移學(xué)習(xí)的方法,即先用模式資料訓(xùn)練神經(jīng)網(wǎng)絡(luò),在此基礎(chǔ)上再用觀測資料對神經(jīng)網(wǎng)絡(luò)進行微調(diào)。做遷移學(xué)習(xí)時,雖然模式資料不能完全反映真實海洋的變化,但是經(jīng)過模式資料的訓(xùn)練,神經(jīng)網(wǎng)絡(luò)可以先收斂到較接近真實的情況,之后,再用觀測資料訓(xùn)練這個神經(jīng)網(wǎng)絡(luò),神經(jīng)網(wǎng)絡(luò)給出的情況就會更加接近真實結(jié)果。遷移學(xué)習(xí)結(jié)合了模式資料數(shù)量多和觀測資料質(zhì)量高兩個優(yōu)點,克服了過往僅采用觀測資料,數(shù)據(jù)集不夠大的問題。實驗二使用觀測資料時,為了正確評估模型的效果,將觀測資料劃分為訓(xùn)練集(1980-1996年)和測試集(1996-2019年)。其中訓(xùn)練集僅用于訓(xùn)練模型,驗證集在后報實驗中用于評估模型效果。
對照實驗三分為兩組,區(qū)別在于輸入變量不同。其中一組只采用SST作為輸入,另一組同時采用SST和HC作為輸入。
除實驗二因需比較不同資料的情況,沒有完全使用全套資料外,實驗一和實驗三均同時使用了模式資料和觀測資料,并同樣采用了遷移學(xué)習(xí)的方法。最終選取在上述實驗中效果最佳的模型進行討論。
圖1給出了各組實驗下模型的預(yù)報技巧,這里使用預(yù)報值與觀測值之間的Pearson相關(guān)系數(shù)作為后報實驗中模型準(zhǔn)確度的度量[19]。從圖1可以看到在任何條件下,相關(guān)系數(shù)均隨著預(yù)報時效的增加而下降。這一點很好理解:海洋運動表現(xiàn)出混沌的特性,初始的微小誤差會隨時間放大,從而影響長時間的預(yù)報效果。
圖1 不同方法下模型的預(yù)報技巧(誤差線表示其他變量變化時,預(yù)報技巧的變化范圍)
圖1(a)展示了單一模型和集合模型的效果差異??梢钥吹?集合模型的引入對于效果有明顯改進,能提升預(yù)報的相關(guān)系數(shù),而且預(yù)報時效越長,提升越明顯。這證明了采用集合模型確實能較大程度地提升預(yù)報效果。
圖1(b)展示了使用不同資料對模型效果的影響??梢钥吹酵瑫r使用模式資料和觀測資料時,模型的表現(xiàn)優(yōu)于單獨使用一種資料。另外可以看到,對于1~7個月(短中期)的預(yù)報,僅用模式資料比僅用觀測資料的效果好,而對于8個月及以上(長期)的預(yù)報,則是僅用觀測資料的效果更好些。這說明模式資料主要改善的是短中期預(yù)報,而觀測資料改善的則是長期預(yù)報。
圖1(c)展示了僅用SST和同時使用SST和HC時效果的差異。引入HC之后,對結(jié)果有負(fù)面的影響。理論上,HC包含了海洋更深層次的信息,應(yīng)該能夠改進預(yù)報的效果,特別是能夠減輕春季預(yù)報障礙[20]。但對神經(jīng)網(wǎng)絡(luò)來說,增加輸入的數(shù)據(jù)量并不是總能提高輸出的效果。例如,Jonathan等[21]在預(yù)報位勢高度場時,對比過僅使用500 hPa的單層資料和同時使用700 hPa、500 hPa、300 hPa的多層資料的效果,后者并沒有表現(xiàn)出比前者更好的改善,研究指出這可能是由于模型沒有輸入變量的垂直結(jié)構(gòu)信息。另外,雖然HC提供了更多的信息,但在實際觀測中,海洋次表層的數(shù)據(jù)比表層稀缺,以及數(shù)值模式對該層次物理過程的刻畫能力,這些都有可能導(dǎo)致HC的數(shù)據(jù)質(zhì)量不如SST??傊?引入HC后預(yù)報效果下降的原因在之后的研究中值得進一步討論。
對照實驗表明,在使用集合模型同時使用模式和觀測資料,以及僅用SST數(shù)據(jù)時,模型的表現(xiàn)最佳。圖2展示了神經(jīng)網(wǎng)絡(luò)模型與傳統(tǒng)動力模式集合平均和統(tǒng)計模式集合平均的預(yù)報技巧比較。可以看到神經(jīng)網(wǎng)絡(luò)模型的效果優(yōu)于傳統(tǒng)動力模式集合平均和統(tǒng)計模式集合平均,在預(yù)報時效為9個月時,相關(guān)系數(shù)仍可以達到0.7以上。
圖2 神經(jīng)網(wǎng)絡(luò)、動力模式集合平均和統(tǒng)計模式集合平均的預(yù)報技巧
圖3比較了在測試時間段(1996-2019年),實際觀測的Nino 3.4指數(shù)和模型給出的預(yù)測值之間的差異。對于短時效的預(yù)報,模型能準(zhǔn)確給出ENSO發(fā)生的時間和強度。隨著預(yù)報時效增加,相關(guān)系數(shù)逐漸下降,集合模型成員的發(fā)散程度增加,對ENSO事件強度的預(yù)報偏小。同時,對長時效的預(yù)報,預(yù)報的峰值的時間存在滯后現(xiàn)象。這種滯后現(xiàn)象被稱為“目標(biāo)時段滯后”[22],在動力模式和統(tǒng)計模式中都存在,而由于包括本文的模型在內(nèi)的統(tǒng)計模式常使用近3個月平均作為輸入,因此與動力模式相比滯后的情況更加突出[22]。另外,由于采用了集合模型,可以給出成員發(fā)散的范圍,這在實際業(yè)務(wù)中具有更多的參考意義。
圖3 不同預(yù)報時效下,觀測的Nino 3.4指數(shù)、模型集合平均預(yù)報和集合成員發(fā)散范圍
圖4給出了在不同的目標(biāo)月份和預(yù)報時效下,模型準(zhǔn)確度的變化。除了準(zhǔn)確度隨預(yù)報時效增加而減小之外,還可以發(fā)現(xiàn)準(zhǔn)確度存在著周年變化。具體而言,在冬春季節(jié)預(yù)測夏季的情況時,準(zhǔn)確率相比于一年中其他時間有明顯的下降,這個ENSO預(yù)報問題被稱為“春季預(yù)報障礙”。進一步計算表明,在預(yù)報時效為6個月時,集合模型50個成員間的方差在春季最大。這說明集合預(yù)報在春季發(fā)散最明顯、預(yù)報的不確定性在春季最大,從另一個角度體現(xiàn)了“春季預(yù)報障礙”。
圖4 不同的目標(biāo)月份和預(yù)報時效對應(yīng)的模型的預(yù)報技巧
春季預(yù)報障礙的存在一般認(rèn)為和季節(jié)循環(huán)、ENSO自身性質(zhì)、初始誤差3個原因有關(guān)[23-24]。南方濤動指數(shù)(SOI)和Nino 3.4指數(shù)的相關(guān)系數(shù)可以表征海-氣相互作用強度。計算了不同月份的SOI和Nino 3.4指數(shù)的相關(guān)系數(shù),發(fā)現(xiàn)相關(guān)性在春季最低,表明海-氣相互作用在春季最不穩(wěn)定,誤差易在春季發(fā)展,這從季節(jié)循環(huán)的角度對春季預(yù)報障礙做出了解釋。
除了模型預(yù)報技巧的年內(nèi)變化外,本文還研究了預(yù)報技巧在年際尺度上的變化(圖5)。統(tǒng)計了在預(yù)報時效為1~6個月時,模型預(yù)報的平均絕對誤差(MAE),計算結(jié)果顯示1~6個月MAE呈現(xiàn)出明顯且一致的年際變化。MAE的年際變化指示出ENSO可預(yù)報性的年際變化,前人的研究指出這與海-氣相互作用的Bjerknes反饋的強度有關(guān)[25]。
根據(jù)國家氣候中心給出的《ENSO歷史事件統(tǒng)計表》,將中等及以上強度的厄爾尼諾發(fā)生的時間段時間繪制在圖上(圖5粉紅色區(qū)域),可以看到在中等及以上強度的厄爾尼諾發(fā)生之前,誤差會有比較大的增長。結(jié)合圖3,可以得出這主要是模型在預(yù)測時對厄爾尼諾的峰值估計偏小、時間估計滯后造成的。另外在對預(yù)報拉尼娜時,誤差增大的現(xiàn)象則表現(xiàn)得不明顯(圖5藍色區(qū)域),呈現(xiàn)出預(yù)報誤差的不對稱性。根據(jù)張雅樂等[26]的研究,在ENSO事件冷暖相位,預(yù)報誤差的增長位相具有顯著的非對稱性,發(fā)展期暖位相預(yù)報誤差強于冷位相,這解釋了預(yù)報誤差不對稱性的現(xiàn)象。
圖5 預(yù)報時效為1~6個月時,模型預(yù)測的平均絕對誤差隨年份的變化
將神經(jīng)網(wǎng)絡(luò)應(yīng)用到ENSO預(yù)報中,利用神經(jīng)網(wǎng)絡(luò)預(yù)測Nino 3.4區(qū)SSTA,并對模型的預(yù)報結(jié)果進行了討論。
對照實驗的結(jié)果表明,采用集合模型的方法不僅能夠提高預(yù)報的準(zhǔn)確度,集合成員給出的發(fā)散范圍在實際業(yè)務(wù)中也具有參考價值。數(shù)值模式的加入,可以彌補實際觀測資料不足的問題。把模式資料和觀測資料相結(jié)合,可以提高模型的預(yù)報技巧。對于預(yù)報因子的選擇,需要通過實驗結(jié)果來判斷其效果。增加預(yù)報因子、擴大輸入量不一定能改善預(yù)報效果。
文中選擇了在多組實驗中表現(xiàn)最佳的模型。該模型在時效為1~9個月的后報實驗中表現(xiàn)優(yōu)于傳統(tǒng)的統(tǒng)計模式集合平均和動力模式集合平均。對后報實驗結(jié)果的分析表明,神經(jīng)網(wǎng)絡(luò)的預(yù)報準(zhǔn)確度存在著年內(nèi)和年際變化。其中,年內(nèi)變化和“春季預(yù)報障礙”有關(guān),年際變化和ENSO可預(yù)報性的年際變化有關(guān)。
神經(jīng)網(wǎng)絡(luò)的預(yù)報也存在著一些不足,例如在做長時效預(yù)報時,對于ENSO事件峰值存在估計偏小、峰值時間預(yù)測滯后的問題。機器學(xué)習(xí)方法在許多領(lǐng)域的應(yīng)用都取得了良好的效果,如何將機器學(xué)習(xí)和氣象有機結(jié)合值得繼續(xù)探索。盡管有時機器學(xué)習(xí)能給出相當(dāng)好的結(jié)果,但與傳統(tǒng)方法相比,神經(jīng)網(wǎng)絡(luò)是一個只負(fù)責(zé)輸入和輸出的“黑箱子”模型。如果能進一步理解神經(jīng)網(wǎng)絡(luò)給出這些結(jié)果的原因,對于更深地認(rèn)識其中的物理機制也是有幫助的。