許柏寧，姜金榮，郝卉群,2，林鵬飛，何丹丹,2

1. 中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心，北京 100190

2. 中國(guó)科學(xué)院大學(xué)，北京 100049

3. 中國(guó)科學(xué)院大氣物理研究所，北京 100029

引言

短期氣候變化包括月尺度、季節(jié)尺度和年際時(shí)間尺度的氣候變化，因其對(duì)人們的生產(chǎn)生活具有重要的影響，而備受人們關(guān)注[1]。其中，年際時(shí)間尺度上的變化又尤其受到人們的關(guān)注[1]。而年際時(shí)間尺度上最強(qiáng)的氣候信號(hào)——厄爾尼諾-南方濤動(dòng) (El Ni?o-Southern Oscillation，ENSO)，自然而然的受到各國(guó)氣象學(xué)家、科學(xué)家的廣泛關(guān)注[1-4]。雖然 ENSO 現(xiàn)象發(fā)生在熱帶太平洋地區(qū)，但是經(jīng)過(guò)海氣的相互作用，其會(huì)以遙相關(guān)的方式間接地對(duì)全球的短期氣候 (氣溫、降水等) 造成重大的影響，引發(fā)多種氣象災(zāi)害，進(jìn)而對(duì)糧食產(chǎn)量造成嚴(yán)重的影響[1-2,4-8]。

鑒于 ENSO 對(duì)全球氣候的重大影響，各國(guó)氣象學(xué)家對(duì) ENSO 展開(kāi)廣泛的研究，以期可以對(duì)其進(jìn)行預(yù)測(cè)。因?yàn)?ENSO 發(fā)生時(shí)的表現(xiàn)為熱帶太平洋區(qū)域的海表面溫度 (Sea Surface Temperature，SST) 的異常，所以，預(yù)報(bào) ENSO 可以等價(jià)于預(yù)測(cè)熱帶太平洋的海表面溫度或海表溫度異常 (Sea Surface Temperature Anomaly，SSTA)。自 20 世紀(jì) 80 年代起，經(jīng)過(guò)三十多年的研究，ENSO 的預(yù)報(bào)模型已經(jīng)可以達(dá)到 6-12 個(gè)月的可信預(yù)報(bào)程度[1]。目前，ENSO 的預(yù)測(cè)模型主要分為兩類：統(tǒng)計(jì)學(xué)模型 (Statistical Model) 和動(dòng)力學(xué)模型 (Dynamical Model)。圖 1 是目前國(guó)際上 14 個(gè)動(dòng)力預(yù)報(bào)模型和 8 個(gè)統(tǒng)計(jì)預(yù)報(bào)模型，自2016年3月以來(lái)對(duì)Ni?o3.4 區(qū)域 (5oS～5oN, 120oW～170oW) 平均海表溫度異常的預(yù)報(bào)結(jié)果。可以看出，雖然模式之間差異比較大，但是模式都能大體上預(yù)測(cè)出 SSTA 的變化趨勢(shì)。此外，最近也有一些學(xué)者考慮到初值的不確定性和模式中物理過(guò)程以及參數(shù)化方案的不確定性，進(jìn)行了 ENSO 集合預(yù)報(bào)的研究[1,9-11]。研究發(fā)現(xiàn)，集合預(yù)報(bào)在預(yù)報(bào)準(zhǔn)確性上要優(yōu)于單一確定性預(yù)報(bào)。

近期，隨著 AI 浪潮的興起，機(jī)器學(xué)習(xí)和深度學(xué)習(xí)又一次進(jìn)入了人們的視線。一些學(xué)者[12-21]開(kāi)始嘗試使用機(jī)器學(xué)習(xí)或深度學(xué)習(xí)的技術(shù)進(jìn)行氣象要素(風(fēng)速、溫度等) 或氣象現(xiàn)象的預(yù)測(cè)，取得了不錯(cuò)的成果，其中，對(duì)于 SST 的預(yù)測(cè)模型[19,20]就可以用來(lái)預(yù)報(bào) ENSO。此外，也有一些學(xué)者[22,23]直接使用機(jī)器學(xué)習(xí)的方法，直接對(duì) ENSO 指標(biāo)進(jìn)行預(yù)測(cè)，也取得了優(yōu)于動(dòng)力預(yù)報(bào)和其他統(tǒng)計(jì)預(yù)報(bào)模式的結(jié)果。但是，目前使用機(jī)器學(xué)習(xí)或深度學(xué)習(xí)對(duì) SST 預(yù)測(cè)的模型中，都是針對(duì)于某個(gè)經(jīng)緯點(diǎn)的時(shí)間序列進(jìn)行預(yù)測(cè)，然后對(duì)某個(gè)區(qū)域的一定網(wǎng)格尺度下所有的經(jīng)緯格點(diǎn)分別運(yùn)行預(yù)測(cè)，進(jìn)而得到整個(gè)區(qū)域的 SST。這樣的做法雖然最后也是得到了某個(gè)區(qū)域的 SST 分布，但是考慮到溫度其實(shí)是一個(gè)同時(shí)具有時(shí)間和空間信息的量，這種單點(diǎn)預(yù)測(cè)的做法明顯忽略了溫度的空間信息，這樣就可能會(huì)破壞掉溫度分布的連續(xù)性。而直接預(yù)報(bào) ENSO 指數(shù)的模型，雖然預(yù)報(bào) ENSO 的結(jié)果要優(yōu)于傳統(tǒng)模式，但是模型最終只是得到了 ENSO 指數(shù)，并沒(méi)有 SST 或SSTA 等氣象要素的信息，這樣的話就無(wú)法根據(jù) SST的分布信息做進(jìn)一步的分析。

本文將以 Ni?o3.4 區(qū)域 ( 5°N～5°S，170°W～120°W) 為研究區(qū)域，用多個(gè)深度學(xué)習(xí)模型對(duì)整個(gè)研究區(qū)域進(jìn)行提前 12 個(gè)月的月均 SSTA 預(yù)測(cè)，根據(jù)月均 SSTA 的預(yù)測(cè)結(jié)果找出最優(yōu)的模型。然后，在最優(yōu)模型的基礎(chǔ)上進(jìn)行一些調(diào)參實(shí)驗(yàn)，找到預(yù)測(cè)效果最優(yōu)的一組參數(shù)。最后，使用在最優(yōu)參數(shù)下訓(xùn)練好的模型，對(duì) Ni?o3.4 區(qū)域進(jìn)行 ENSO 指數(shù) (Ni?o3.4 指數(shù))的預(yù)測(cè)。

圖1 國(guó)際上 14 個(gè)動(dòng)力預(yù)報(bào)模式和 8 個(gè)統(tǒng)計(jì)預(yù)報(bào)模式的預(yù)報(bào)結(jié)果 (圖片來(lái)源：https://iri.columbia.edu/our-expertise/climate/forecasts/enso/current/?enso_tab=enso-sst_table)Fig. 1 The prediction results of international fourteen dynamical models and eight statistical models (Figure source: https://iri.columbia.edu/our-expertise/climate/forecasts/enso/current/?enso_tab=enso-sst_table)

1 ENSO 指數(shù)和模型關(guān)鍵層簡(jiǎn)介

1.1 ENSO 指數(shù)

ENSO 是一個(gè)發(fā)生在熱帶太平洋并且涉及到海氣不同方面的復(fù)雜過(guò)程，因此，有很多在不同方面描述 ENSO 事件的指數(shù)：(1) 從氣壓方面衡量的南方濤動(dòng)指數(shù)[24](Southern Oscillation Index，簡(jiǎn)稱 SOI)以及改進(jìn)的赤道南方濤動(dòng)指數(shù) (Equatorial Southern Oscillation Index，簡(jiǎn)稱 ESOI)；(2) 根據(jù)不同海域的 SST 方面進(jìn)行衡量的 Ni?o 指數(shù)，包括 Ni?o1+2、Ni?o3、Ni?o4 指數(shù)，以及還有被認(rèn)為更能衡量ENSO 的 Ni?o3.4 指數(shù)[25]；(3) 根據(jù)長(zhǎng)波輻射進(jìn)行衡量 ENSO 的指標(biāo)[26]； (4) 利用風(fēng)進(jìn)行衡量 ENSO 的指標(biāo)。本實(shí)驗(yàn)使用常用的 Ni?o3.4 作為衡量 ENSO 強(qiáng)度的指標(biāo)，而 Ni?o3.4 指標(biāo)是通過(guò)計(jì)算 Ni?o3.4 區(qū)域的SSTA 平均值得到的。

1.2 模型網(wǎng)絡(luò)層簡(jiǎn)介

本實(shí)驗(yàn)所使用的神經(jīng)網(wǎng)絡(luò)層主要包括卷積神經(jīng)網(wǎng)絡(luò)[27-28](Convolutional Neural Network，CNN) 中的卷積層和卷積長(zhǎng)短時(shí)記憶[29](Convolutional Long Short-Term Memory，ConvLSTM) 單元。下面，對(duì)兩者進(jìn)行簡(jiǎn)單的介紹。

1.2.1 卷積神經(jīng)網(wǎng)絡(luò)

卷積神經(jīng)網(wǎng)絡(luò)的提出與動(dòng)物體的視覺(jué)處理有密切的關(guān)系。Hubel 等[30]學(xué)者通過(guò)研究貓和猴子的視覺(jué)神經(jīng)后，提出了視覺(jué)感受野的概念。K. Fukushima等[31]學(xué)者根據(jù)視覺(jué)感受野的概念提出的神經(jīng)認(rèn)知機(jī) (Neocognitron) 可以認(rèn)為是 CNN 的雛形。后來(lái)，LeCun 等[27,32]學(xué)者借用神經(jīng)認(rèn)知機(jī)的思路，利用后向傳播 (Backpropagation，簡(jiǎn)稱 BP) 算法，設(shè)計(jì)并訓(xùn)練了 CNN。后來(lái)，Krizhevsky 等[33]實(shí)現(xiàn)的 CNN 模型在ImageNet 的圖像識(shí)別比賽中大放異彩后，CNN 開(kāi)始被人們熟知，并大量應(yīng)用于圖像分割、圖像分類等圖像類任務(wù)中。后來(lái)，CNN 因其卷積過(guò)程類似于自然語(yǔ)言處理中的 n-gram，也被一些學(xué)者[34-36]逐漸應(yīng)用到自然語(yǔ)言處理領(lǐng)域中。

CNN 中一般都包含有卷積層和下采樣 (pooling)層。其中，卷積層中完成卷積運(yùn)算，而 pooling 層主要用來(lái)保證不變性 (平移不變性、旋轉(zhuǎn)不變性等) 以及減少輸出大小從而減少計(jì)算量。2 維卷積運(yùn)算的前兩次操作如圖 2 所示，其中 A 為 5×5 的二維輸入，W為 3×3 的卷積核 (kernel)，B 為卷積后的輸出 (3×3)，一般稱為特征映射 (feature map)，兩個(gè)方向 (橫向和縱向) 上的卷積步長(zhǎng)均為 1。具體過(guò)程為：(1) 用 W 對(duì) A的左上角 3×3 的矩陣 (淺橙色) 進(jìn)行按位相乘并求和的運(yùn)算后，得到 B11，完成第一次操作，如圖 2 上半部所示；(2) 向右移動(dòng) 1 個(gè)單位 (即橫向的步長(zhǎng)大小)，用 W 對(duì) A 中上部的 3×3 的矩陣 (淺黃色) 進(jìn)行按位相乘并求和的運(yùn)算后，得到 B12，完成第二次操作，如圖 2 下半部所示；以此類推完成整個(gè)卷積運(yùn)算，直至得到 B。

1.2.2 ConvLSTM 單元 [29]

正如 ConvLSTM 的名字所示，ConvLSTM 可以被定義成一種結(jié)合了卷積的長(zhǎng)短時(shí)記憶 (Long Short-Term Memory，簡(jiǎn)稱 LSTM) 單元 (cell)。近年來(lái)，以 LSTM 單元或 GRU 單元為主體的循環(huán)神經(jīng)網(wǎng)絡(luò)(Recurrent Neural Network，簡(jiǎn)稱 RNN) 在自然語(yǔ)言處理等需要處理時(shí)序相關(guān)問(wèn)題的領(lǐng)域中取得很大的成功。但是，不管是傳統(tǒng)簡(jiǎn)單的 RNN 中的單元，還是改進(jìn)過(guò)后的 LSTM 單元、GRU 單元，它們?cè)诿恳粋€(gè)時(shí)間步可接收的輸入都是向量或標(biāo)量，無(wú)法接收矩陣以及維度更高的數(shù)據(jù)。而 ConvLSTM 的提出，使得RNN 網(wǎng)絡(luò)可以在每個(gè)時(shí)間步接收矩陣以及更高維度的數(shù)據(jù)作為輸入。這是因?yàn)?ConvLSTM 把傳統(tǒng) LSTM中一些矩陣乘向量的運(yùn)算換成了卷積核對(duì)輸入數(shù)據(jù)的卷積運(yùn)算，而這樣的話，因?yàn)榫矸e操作可以適用于任意維度的矩陣，所以 ConvLSTM 單元的每個(gè)時(shí)間步就可以接收維度大于 1 的數(shù)據(jù)，進(jìn)而可以對(duì)輸入數(shù)據(jù)同時(shí)進(jìn)行時(shí)空兩個(gè)維度上的特征學(xué)習(xí)。ConvLSTM 相對(duì)于 LSTM 的具體改變?nèi)鐖D 3 中公式所示，其中，x和h為向量，X和H則不限維度，* 為哈達(dá)馬積(Hadamard product)， 為卷積運(yùn)算。

圖3 LSTM 到 ConvLSTM 公式的改變 (左：LSTM；右：ConvLSTM)Fig. 3 Formula changes from LSTM to ConvLSTM (left:LSTM; right: ConvLSTM)

2 數(shù)據(jù)

本實(shí)驗(yàn)中所使用的數(shù)據(jù)主要包括 SSTA、表面壓力 (P_SURF)、經(jīng)向風(fēng)速 (V_WIND)、緯向風(fēng)速 (U_WIND)、經(jīng)向流速 (V_CURRENT) 和緯向流速 (U_CURRENT)。其中，SSTA 來(lái)源于 NOAA OI SST V2 High Resolution Dataset[37]數(shù)據(jù)集，剩下的 5 種變量均來(lái)源于 NCEP Climate Forecast System Reanalysis (CFSR) 6-hourly Products[38,39]和 NCEP Climate Forecast System Version 2 (CFSv2) 6-hourly Products[40]數(shù)據(jù)集。因?yàn)閿?shù)據(jù)來(lái)源不同，導(dǎo)致網(wǎng)格的分辨率以及時(shí)間范圍均有一些差異 (例如 SSTA是 0.25°*0.25° 網(wǎng)格的數(shù)據(jù)，而 P 和 U_WIND 均為0.5°*0.5° 網(wǎng)格的數(shù)據(jù))，對(duì)數(shù)據(jù)進(jìn)行篩選和預(yù)處理后，得到了 Ni?o3.4 區(qū)域上 1981 年 9 月 1 日至 2017年 11 月 25 日共 13235 個(gè)日均數(shù)據(jù)，數(shù)據(jù)的分辨率為 0.5°*0.5°。

考慮到 ENSO 的預(yù)測(cè)一般都是基于周均 SSTA或月均 SSTA 的預(yù)測(cè)，所以本實(shí)驗(yàn)擬進(jìn)行月均SSTA 的預(yù)測(cè)，進(jìn)而計(jì)算出 Ni?o3.4 指數(shù)，預(yù)測(cè)ENSO 現(xiàn)象。而通過(guò)現(xiàn)有的日均數(shù)據(jù)只能得到 535個(gè)真實(shí)月均數(shù)據(jù) (1981 年 9 月 1 日至 2017 年 11 月共 535 個(gè)月)，但是，這個(gè)數(shù)據(jù)規(guī)模顯然是不足以支撐深度學(xué)習(xí)模型的訓(xùn)練。因此，本文采用 29 天滑動(dòng)平均的方法進(jìn)行月均數(shù)據(jù)的擴(kuò)增。具體做法為：在所有日均數(shù)據(jù) (d1，d2，d3，…，d13235) 上，滑動(dòng)一個(gè)時(shí)間長(zhǎng)度為 29 的窗口，對(duì)每個(gè)窗口內(nèi)的數(shù)據(jù)計(jì)算平均值作為窗口正中間的數(shù)據(jù)所在月份的近似月均數(shù)據(jù)。例如，第一個(gè)滑動(dòng)窗口中的數(shù)據(jù)為 d1，d2，d3，…，d29，對(duì)這 29 個(gè)數(shù)據(jù)求均值后，得到的數(shù)據(jù)作為 d15所在月份的近似月均數(shù)據(jù)，以此類推可以得到 d16，d17，…，d13221所在月份的近似月均數(shù)據(jù)。經(jīng)過(guò)數(shù)據(jù)集擴(kuò)增后，總共可以得到 13207 個(gè)近似月均數(shù)據(jù)。經(jīng)過(guò)分析，擴(kuò)增月均數(shù)據(jù)中超過(guò) 60%的數(shù)據(jù)點(diǎn)與真實(shí)月均數(shù)據(jù)的絕對(duì)誤差小于 0.1℃，超過(guò) 82% 的數(shù)據(jù)點(diǎn)的絕對(duì)誤差小于 0.2℃，約 98% 的數(shù)據(jù)點(diǎn)的絕對(duì)誤差小于 0.5℃，因此擴(kuò)增數(shù)據(jù)具有有效性。

3 模型結(jié)構(gòu)設(shè)計(jì)及實(shí)驗(yàn)過(guò)程

3.1 模型結(jié)構(gòu)

本實(shí)驗(yàn)采用兩類訓(xùn)練模型——Sequence-to-One(Seq2One) 和 Sequence-to-Sequence (Seq2Seq) 進(jìn)行訓(xùn)練。兩類訓(xùn)練模型主要有兩個(gè)不同之處：

(1) 輸出長(zhǎng)度不同。兩種訓(xùn)練模型的輸入均為一定時(shí)間長(zhǎng)度 (如 12 個(gè)月) 的時(shí)空序列，而輸出分別為單一時(shí)間點(diǎn)的 SSTA 數(shù)據(jù) (即序列的時(shí)間長(zhǎng)度為 1) 和一定時(shí)間長(zhǎng)度 (如 12 個(gè)月) SSTA 時(shí)空序列。例如，輸入均為 2011 年 1 月至 12 月長(zhǎng)度為 12 個(gè)月的時(shí)空序列，對(duì)于 Sequence-to-One 模型，輸出可以為 2012年 1 月到 12 月某一個(gè)月份的月均 SSTA 數(shù)據(jù)，而對(duì)于 Sequence-to-Sequence 模型，輸出可以為 2012 年 1月到 12 月的月均 SSTA 時(shí)空序列。

(2) 進(jìn)行預(yù)測(cè)的形式不同。對(duì)于 Sequence-to-One模型，不同的時(shí)間間隔 (要預(yù)測(cè)的月份與輸入序列的最后一個(gè)月份間的差) 對(duì)應(yīng)于不同的模型，如在上述例子中，對(duì) 2012 年 1 月到 12 月的預(yù)測(cè)就分別對(duì)應(yīng)于12 個(gè)模型；而對(duì)于 Sequence-to-Sequence 模型，則是根據(jù)起報(bào)月 (輸出序列的起始月份) 不同而訓(xùn)練不同的模型進(jìn)行預(yù)測(cè)。

本實(shí)驗(yàn)中設(shè)計(jì)使用了 3 種網(wǎng)絡(luò)結(jié)構(gòu)不同的Sequence-to-One 模型和 1 種 Sequence-to-Sequence模型。具體的網(wǎng)絡(luò)結(jié)構(gòu)示意圖如圖 4 和圖 5 所示。其中@前為 kernel 數(shù)目，@ 后為 kernel 的大小，Sequence-to-Sequence 模型中的 ? 表示 kernel 的個(gè)數(shù)會(huì)隨著輸入數(shù)據(jù)的變量個(gè)數(shù) (channel) 的變化而變化。

不同的網(wǎng)絡(luò)結(jié)構(gòu)中，一般來(lái)說(shuō)由 CNN 進(jìn)行空間特征的學(xué)習(xí)，ConvLSTM 進(jìn)行時(shí)間維度的學(xué)習(xí)，考慮到 ConvLSTM 中也包括卷積過(guò)程以及 3維 CNN 也可以同時(shí)進(jìn)行時(shí)間維度上的卷積，所以也有只使用 ConvLSTM (圖 4.b) 和 CNN (圖 4.a) 的網(wǎng)絡(luò)結(jié)構(gòu)。Sequence-to-Sequence 模型的網(wǎng)絡(luò)結(jié)構(gòu)則借用自然語(yǔ)言處理領(lǐng)域中的 Seq2Seq 結(jié)構(gòu)，只是把 Seq2Seq 結(jié)構(gòu)中的 RNN 相關(guān)的單元換成了ConvLSTM。

圖4 Sequence-to-One 模型示意圖 (a.CNNs 模型，b.ConvLSTMs 模型，c. CNNs+ConvLSTMs 模型)Fig. 4 Schematic fi gure of Sequence-to-One Models (a. CNNs model, b. ConvLSTMs model, c. CNNs+ConvLSTMs model)

圖5 Sequence-to-Sequence 模型示意圖Fig. 5 Schematic fi gure of Sequence-to-Sequence model

3.2 模型訓(xùn)練

本文把 SSTA 的預(yù)測(cè)問(wèn)題轉(zhuǎn)換為時(shí)空序列的預(yù)測(cè)問(wèn)題，通過(guò)使用大小為序列長(zhǎng)度的滑動(dòng)窗口，在擴(kuò)增的月均數(shù)據(jù)集上滑動(dòng)出模型訓(xùn)練所需的樣本，得到的樣本根據(jù) 8:2 的比例劃分為訓(xùn)練集和測(cè)試集。此外，本文采用均方根誤差 (Root Mean Square Error，RMSE) 作為訓(xùn)練過(guò)程中的損失函數(shù)。

3.3 實(shí)驗(yàn)過(guò)程

本文的實(shí)驗(yàn)包括五組實(shí)驗(yàn)，其中，第一組實(shí)驗(yàn)用來(lái)確定最優(yōu)的預(yù)測(cè)模型，第二組實(shí)驗(yàn)用來(lái)確定最優(yōu)的一組超參數(shù)，第三組實(shí)驗(yàn)用來(lái)確定最優(yōu)變量個(gè)數(shù)，剩余兩組實(shí)驗(yàn)用來(lái)與其他 ENSO 預(yù)報(bào)模型進(jìn)行對(duì)比。具體的實(shí)驗(yàn)如下：

(1) 確定最優(yōu)模型結(jié)構(gòu)。使用上述設(shè)計(jì)的四種模型結(jié)構(gòu)，在 Ni?o3.4 區(qū)域上進(jìn)行實(shí)驗(yàn)，通過(guò) RMSE 的對(duì)比找出最優(yōu)的模型結(jié)構(gòu)；

(2) 確定在最優(yōu)的模型基礎(chǔ)上的一組最優(yōu)的超參數(shù)。調(diào)節(jié)優(yōu)化器種類 (Adam、Adadelta、SGD、Adgrad)、是否引入正則化等超參數(shù)，選擇一組最優(yōu)的超參數(shù)。

(3) 確定最佳的輸入變量的個(gè)數(shù)。分別使用SSTA、SSTA+P+U_WIND (P 和 U_WIND 與 SSTA 的相關(guān)性略強(qiáng)) 和所有變量作為輸入進(jìn)行三次實(shí)驗(yàn)，以確定最佳的輸入變量個(gè)數(shù)；

(4) 用上述實(shí)驗(yàn)確定的最優(yōu)的模型和參數(shù) (優(yōu)化器種類、是否正則化、輸入變量個(gè)數(shù)等) ，訓(xùn)練出用于預(yù)測(cè) SSTA 的模型，與動(dòng)力學(xué)模型[9]進(jìn)行對(duì)比，以檢驗(yàn)?zāi)Ｐ皖A(yù)測(cè)效果的好壞；

(5) 使用訓(xùn)練好的模型進(jìn)行 Ni?o3.4 區(qū)域的 SSTA預(yù)測(cè) (即 Ni?o3.4 指數(shù)預(yù)測(cè)) 實(shí)驗(yàn)，與其他動(dòng)力學(xué)模型和統(tǒng)計(jì)模型的預(yù)測(cè)結(jié)果進(jìn)行對(duì)比。

需要說(shuō)明的是，因?yàn)?SSTA 的序列自相關(guān)性最強(qiáng)，因此，在前兩組實(shí)驗(yàn)中只使用 SSTA 作為模型輸入。

4 實(shí)驗(yàn)結(jié)果與分析

第一組實(shí)驗(yàn)的實(shí)驗(yàn)結(jié)果如圖 6。從圖 6 中不難看出，就 RMSE 的對(duì)比 (圖 6.a) 而言，Seq2Seq 模型不論在訓(xùn)練集還是測(cè)試集上都要優(yōu)于 Seq2One 模型。在相關(guān)系數(shù) (圖 6.b) 的對(duì)比上，中短期 (lead months ＜6) 預(yù)測(cè)中，Seq2One 模型相對(duì)表現(xiàn)更好一點(diǎn)，而中長(zhǎng)期 (lead months ＞ 6) 預(yù)測(cè)中反而是 Seq2Seq 模型表現(xiàn)更好一點(diǎn)。其中，Seq2Seq 模型在向后預(yù)報(bào)到 3-6 個(gè)月 (lead months = 3-6) 時(shí)，相關(guān)系數(shù)有明顯的降低，這與動(dòng)力學(xué)預(yù)報(bào)上的“春季預(yù)報(bào)障礙”現(xiàn)象比較吻合，這可能就說(shuō)明 Seq2Seq 模型的預(yù)測(cè)過(guò)程類似于動(dòng)力學(xué)預(yù)報(bào)。而 Seq2One 模型 (除 CNNs 模型之外) 則沒(méi)有類似“春季預(yù)報(bào)障礙”的現(xiàn)象。但是，整體來(lái)看，Seq2Seq 模型還是要優(yōu)于 Seq2One 模型的，而且ENSO 預(yù)報(bào)一般要越能提前并準(zhǔn)確地預(yù)報(bào)越好，在中長(zhǎng)期的預(yù)報(bào)技巧 (相關(guān)系數(shù)) 上，顯然也是 Seq2Seq 模型更占優(yōu)勢(shì)。因此，綜合 RMSE 和相關(guān)系數(shù)的對(duì)比，可以看出 Seq2Seq 模型是設(shè)計(jì)的四種模型結(jié)構(gòu)中最優(yōu)的。此外，值得注意的是，四種模型在訓(xùn)練集 (實(shí)線)和測(cè)試集 (虛線) 上的表現(xiàn)有一定的差距。

圖6 RMSE 和相關(guān)系數(shù)對(duì)比結(jié)果圖 (train：訓(xùn)練集；test：測(cè)試集)Fig. 6 RMSE and correlation coefficient comparison result fi gure (train: training set; test: test set)

圖7 不同超參數(shù)對(duì)比實(shí)驗(yàn)結(jié)果圖 (測(cè)試集)Fig. 7 Results of different hyperparametric contrast experiments (test set)

圖7為不同超參數(shù)實(shí)驗(yàn)在測(cè)試集上的實(shí)驗(yàn)結(jié)果的對(duì)比圖，其中考慮到 Seq2Seq 模型的編碼器 (Encoder)部分對(duì)輸入序列可以有三種方式——正序輸入、逆序輸入和使用雙向編碼器 (Bi-Encoder) 接收雙向輸入，所以也把編碼器不同的輸入方式也考慮為超參數(shù)之一。首先，就不同的編碼器輸入方式 (圖 7.a 和圖 7.b) 而言，在 RMSE 的對(duì)比上，正序輸入 (normal_encoder_input) 只有在短期 (lead months ＜ 5) 預(yù)測(cè)時(shí)表現(xiàn)得略差，在中長(zhǎng)期 (lead months ≥ 5) 預(yù)測(cè)時(shí)要明顯好于其他兩種輸入方式；而在相關(guān)系數(shù)的對(duì)比上，三種輸入在短期 (lead months ＜ 5) 預(yù)測(cè)時(shí)的表現(xiàn)相差無(wú)幾，而在中長(zhǎng)期 (lead months ≥ 5) 預(yù)測(cè)時(shí)，隨著提前預(yù)報(bào)月份 (lead months) 的增加，正序輸入與其他兩種輸入的差距逐漸增大，最終，正序輸入的相關(guān)系數(shù)比其他兩種輸入的相關(guān)系數(shù)要高 0.1 左右。因此，編碼器使用正序輸入可以帶來(lái)更好的預(yù)測(cè)結(jié)果。其次，就是否使用 L2 正則化 (圖 7.c 和圖 7.d) 這一點(diǎn)來(lái)說(shuō)，不使用 L2 正則化的模型的結(jié)果 (包括 RMSE 和相關(guān)系數(shù)) 都要略好于使用 L2 正則化的模型的結(jié)果。其中，短期 (lead months ＜ 3) 預(yù)測(cè)時(shí)，不使用 L2 正則化的模型的RMSE更是明顯好于使用 L2 正則化的模型的RMSE。因此，不使用L2正則化更有利于模型得到好的預(yù)測(cè)結(jié)果。最后，針對(duì)于不同的優(yōu)化器 (圖 7.e 和圖 7.f)，Adam 優(yōu)化器在 RMSE 和相關(guān)系數(shù)上都要優(yōu)于其他三種優(yōu)化器。綜上所述，對(duì)于 Seq2Seq 模型而言，最優(yōu)的超參數(shù)為：編碼器使用正序輸入、不使用L2 正則化和使用 Adam 優(yōu)化器。

第三組實(shí)驗(yàn)的實(shí)驗(yàn)結(jié)果如圖 8 所示。可以看出增加變量個(gè)數(shù)只能降低訓(xùn)練集上的 RMSE，而對(duì)測(cè)試集來(lái)說(shuō)基本沒(méi)有提升，甚至中長(zhǎng)期 (lead months ＞ 5)的預(yù)測(cè)反而是只使用 SSTA 要更好。因此，這就說(shuō)明只使用SSTA做模型的輸入就可以達(dá)到比較好的預(yù)測(cè)結(jié)果。此外，使用 3 個(gè)變量 (SSTA、P_SURF、U_WIND) 作為輸入和使用 6 個(gè)變量作為輸入對(duì)模型的預(yù)測(cè)結(jié)果基本沒(méi)有影響。這可能說(shuō)明增加的 3 個(gè)變量(V_WIND、U_CURRENT、V_CURRENT) 對(duì) SSTA的預(yù)測(cè)基本沒(méi)有影響。

使用 Seq2Seq 模型，用 SSTA 作為輸入變量，用前述實(shí)驗(yàn)所確定的超參數(shù) (編碼器使用正序輸入、不使用 L2 正則化和使用 Adam 優(yōu)化器) 進(jìn)行 SSTA 的預(yù)測(cè)實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果如圖 9 所示。其中，圖中陰影區(qū)域表示 12 個(gè)模型 (分別以 12 個(gè)月作為起報(bào)月得到的模型) 在測(cè)試集上的結(jié)果的分布范圍，而紅色實(shí)線為 12個(gè)模型結(jié)果的平均值，即不考慮不同起報(bào)月預(yù)測(cè)結(jié)果的差別。由圖 9 可以看出，就 RMSE 來(lái)看，Seq2Seq模型在短期 (lead months ＜ 3) 預(yù)測(cè)時(shí)的結(jié)果要略差于單一動(dòng)力學(xué)模型 (deterministic model) 和集合預(yù)報(bào)(ensemble model) 的結(jié)果。但是隨著預(yù)測(cè)超前的月份 (lead months) 的增加，Seq2Seq 模型逐漸顯示出優(yōu)勢(shì)。當(dāng) lead months ≥ 3 時(shí)，Seq2Seq 模型的預(yù)測(cè)結(jié)果就已經(jīng)優(yōu)于單一動(dòng)力學(xué)模型的預(yù)測(cè)結(jié)果了，而當(dāng) lead months ＞ 7 后，其預(yù)測(cè)結(jié)果更是優(yōu)于集合預(yù)報(bào)的預(yù)測(cè)結(jié)果。但是，就相關(guān)系數(shù)而言，Seq2Seq 模型與動(dòng)力學(xué)模型還是有一定的差距的。此外，經(jīng)過(guò)分析，測(cè)試集上的較差的邊界 (RMSE 中的上邊界和相關(guān)系數(shù)的下邊界) 應(yīng)該與“春季預(yù)報(bào)障礙”現(xiàn)象有一定關(guān)系。

圖8 不同變量個(gè)數(shù)實(shí)驗(yàn)結(jié)果圖Fig. 8 Results of different number of variables’ experiments

圖9 Seq2Seq 模型與動(dòng)力學(xué)模型對(duì)比Fig. 9 Comparison of Seq2Seq model and dynamical models

最后，利用訓(xùn)練好的 Seq2Seq 模型對(duì) 2015 之后Ni?o3.4 區(qū)域的 SSTA 進(jìn)行預(yù)測(cè)實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果如圖10 所示。從圖中可以看出，Seq2Seq 模型對(duì) 2015 年出現(xiàn)的強(qiáng)厄爾尼諾現(xiàn)象擬合的較差。但是，相較于其他 ENSO 預(yù)報(bào)模型的預(yù)測(cè)結(jié)果 (圖 10 虛線，數(shù)據(jù)來(lái)源于 https://iri.columbia.edu/our-expertise/climate/enso/)，Seq2Seq 模型的預(yù)測(cè)結(jié)果還算可以，而且在 2017 年 5月左右的峰后可以較好地預(yù)測(cè)出回落。因此，整體上來(lái)看，目前的 Seq2Seq 模型可以較好地預(yù)測(cè)出 SSTA的趨勢(shì)，但是峰值強(qiáng)度的預(yù)測(cè)還有一定的差距。

5 結(jié)語(yǔ)

本文通過(guò)一組對(duì)比實(shí)驗(yàn)，從四種設(shè)計(jì)的網(wǎng)絡(luò)模型中挑選了一種可以用來(lái) SSTA 預(yù)測(cè)的模型——Seq2Seq 模型。隨后又確定了可以讓 Seq2Seq 模型取得最優(yōu)預(yù)測(cè)結(jié)果的一組超參數(shù)和輸入變量個(gè)數(shù)。然后，用此模型與傳統(tǒng)的動(dòng)力學(xué)模型進(jìn)行了 RMSE 和預(yù)報(bào)技巧 (相關(guān)系數(shù)) 的對(duì)比，結(jié)果表明 Seq2Seq 模型在中長(zhǎng)期 (lead months ＞ 7) 的預(yù)測(cè)結(jié)果要優(yōu)于動(dòng)力學(xué)模型，而在預(yù)報(bào)技巧上的表現(xiàn)要略差于動(dòng)力學(xué)模型，但是中長(zhǎng)期較好的預(yù)測(cè)結(jié)果更有利于 Seq2Seq 模型用于 ENSO 現(xiàn)象的預(yù)報(bào)。最后，使用訓(xùn)練好的Seq2Seq模型對(duì) 2015 年 1 月到 2018 年 8 月的 SSTA 進(jìn)行了預(yù)測(cè)，結(jié)果表明 Seq2Seq 模型可以較好地預(yù)測(cè)出 SSTA的變化趨勢(shì)，但是預(yù)測(cè)的峰值強(qiáng)度 (特別是厄爾尼諾的強(qiáng)度) 要弱于觀測(cè)值。此外，與其他ENSO預(yù)報(bào)模型相比，Seq2Seq 模型的表現(xiàn)處于中等水平。

當(dāng)然，目前的 Seq2Seq 模型也存在一些問(wèn)題：(1)測(cè)試集上有較為明顯的類似于動(dòng)力學(xué)模型中的“春季預(yù)報(bào)障礙”的現(xiàn)象；(2) 訓(xùn)練集與測(cè)試集有較大差距；(3) 峰值強(qiáng)度較低。

因此，在今后的實(shí)驗(yàn)中可以嘗試解決上述問(wèn)題，以優(yōu)化現(xiàn)有的模型。本文中只是進(jìn)行了相對(duì)簡(jiǎn)單的調(diào)參工作，沒(méi)有對(duì)正則化的懲罰因子的大小，不同優(yōu)化器的學(xué)習(xí)率等進(jìn)行更深入的調(diào)參工作，隨著今后調(diào)參工作的進(jìn)一步深入，或許可以找到能減輕訓(xùn)練集中類 “春季預(yù)報(bào)障礙”的現(xiàn)象的超參數(shù)以及縮小訓(xùn)練集和測(cè)試集之間的差距。而且，也可以嘗試?yán)肧eq2One 模型，特別是 CNNs_ConvLSTMs 模型，減輕 Seq2Seq 模型在訓(xùn)練集中的類 “春季預(yù)報(bào)障礙”現(xiàn)象。針對(duì)峰值強(qiáng)度較低的問(wèn)題，或許可以嘗試引入一些隨機(jī)擾動(dòng)的機(jī)制讓模型可以預(yù)測(cè)出更大的強(qiáng)度。此外，因?yàn)楸疚闹械臄?shù)據(jù)量比較小，導(dǎo)致無(wú)法支撐較大模型的訓(xùn)練，如果可以得到更充足的數(shù)據(jù)，或許可以訓(xùn)練出預(yù)測(cè)效果更好的模型。

圖10 預(yù)測(cè)實(shí)驗(yàn)結(jié)果對(duì)比圖 (虛線為其他 ENSO 預(yù)報(bào)模型的預(yù)測(cè)結(jié)果)Fig. 10 Comparison of prediction experiments’ results (dotted line is the prediction result of other ENSO prediction models)

致謝

本研究受到國(guó)家重點(diǎn)研發(fā)計(jì)劃 (2016YFB0200800)，重點(diǎn)自然科學(xué)基金重點(diǎn)基金 (91530324) 和中國(guó)科學(xué)院知識(shí)創(chuàng)新項(xiàng)目 (XXH13504-03-02) 的資助，在此特別感謝。

[1]鄭飛. ENSO 集合預(yù)報(bào)研究 [D]; 中國(guó)科學(xué)院研究生院(大氣物理研究所), 2007.

[2]鄭冬曉, 楊曉光. ENSO 對(duì)全球及中國(guó)農(nóng)業(yè)氣象災(zāi)害和糧食產(chǎn)量影響研究進(jìn)展 [J]. 氣象與環(huán)境科學(xué), 2014,37(4): 90-101.

[3]李曉燕, 翟盤茂. ENSO 事件指數(shù)與指標(biāo)研究 [J]. 氣象學(xué)報(bào), 2000, 58(1): 102-109.

[4]任福民, 袁媛, 孫丞虎, 曹璐. 近 30 年 ENSO 研究進(jìn)展回顧 [J]. 氣象科技進(jìn)展, 2012, 02(3): 17-24.

[5]許武成, 馬勁松, 王文. 關(guān)于 ENSO 事件及其對(duì)中國(guó)氣候影響研究的綜述 [J]. 氣象科學(xué), 2005, 25(2): 212-220.

[6]秦堅(jiān)肇. 兩類厄爾尼諾對(duì)次年夏季中國(guó)東部降水的可能影響 [D]; 中國(guó)氣象科學(xué)研究院, 2014.

[7]吳萍, 丁一匯, 柳艷菊. 厄爾尼諾事件對(duì)中國(guó)夏季水汽輸送和降水分布影響的新研究[J]. 氣象學(xué)報(bào), 2017, 75(3):371-383.

[8]袁媛, 楊輝, 李崇銀. 不同分布型厄爾尼諾事件及對(duì)中國(guó)次年夏季降水的可能影響 [J]. 氣象學(xué)報(bào), 2012, 70(3):467-478.

[9]ZHENG F, ZHU J. Improved ensemble-mean forecasting of ENSO events by a zero-mean stochastic error model of an intermediate coupled model[J]. Climate Dynamics,2016, 47(12): 1-15. 10.1007/s00382-016-3048-0.

[10]王輝, 劉娜, 李本霞, 李響. 海洋可預(yù)報(bào)性和集合預(yù)報(bào)研究綜述 [J]. 地球科學(xué)進(jìn)展, 2014, 29(11): 1212-1225.

[11]周廣慶, 李旭, 曾慶存. 一個(gè)可供 ENSO 預(yù)測(cè)的海氣耦合環(huán)流模式及 1997/1998ENSO 的預(yù)測(cè) [J]. 氣候與環(huán)境研究, 1998, 29(4): 62-70.

[12]SALMAN A G, KANIGORO B, HERYADI Y. Weather forecasting using deep learning techniques; proceedings of the Advanced Computer Science and Information Systems(ICACSIS), 2015 International Conference on, F, 2015[C]. IEEE. 281-285.

[13]SOMAN S S, ZAREIPOUR H, MALIK O, MANDAL P. A review of wind power and wind speed forecasting methods with different time horizons; proceedings of the North American Power Symposium, F, 2010 [C]. 1-8.

[14]劉鑫達(dá). 基于深度學(xué)習(xí)的氣象溫度預(yù)測(cè)研究 [D]; 寧夏大學(xué), 2016.

[15]燕東渭, 孫田文, 楊艷, 方建剛, 劉志鏡. 支持向量數(shù)據(jù)描述在西北暴雨預(yù)報(bào)中的應(yīng)用試驗(yàn) [J]. 應(yīng)用氣象學(xué)報(bào),2007, 18(5): 676-681.

[16]張振. 基于深度學(xué)習(xí)的全天空云圖分類方法研究 [D];北京交通大學(xué), 2017.

[17]朱磊. 基于圖像的結(jié)露與結(jié)霜現(xiàn)象自動(dòng)化檢測(cè)技術(shù)研究[D]; 華中科技大學(xué), 2015.

[18]FENG Q Y, VASILE R, SEGOND M, GOZOLCHIANI A, WANG Y, ABEL M, HAVLIN S, BUNDE A,DIJKSTRA H A. ClimateLearn: A machine-learning approach for climate prediction using network measures[J].Geoscientific Model Development Discussions, 2016:1-18.

[19]ZHANG Q, WANG H, DONG J Y, ZHONG G Q, SUN X.Prediction of Sea Surface Temperature Using Long Short-Term Memory[J]. Ieee Geosci Remote S, 2017, 14(10):1745-1749. 10.1109/Lgrs.2017.2733548.

[20]AGUILARMARTINEZ S, HSIEH W W. Forecasts of Tropical Pacific Sea Surface Temperatures by Neural Networks and Support Vector Regression[J]. International Journal of Oceanography, 2014, 2009(1687-9406)

[21]KALRA A, AHMAD S. Using oceanic-atmospheric oscillations for long lead time stream fl ow forecasting[J].Water Resources Research, 2009, 45(3): 450-455.

[22]LIMA C H R, LALL U, JEBARA T, BARNSTON A G. Machine Learning Methods for ENSO Analysis and Prediction [M]//LAKSHMANAN V, GILLELAND E, MCGOVERN A, et al. Machine Learning and Data Mining Approaches to Climate Science: Proceedings of the 4th International Workshop on Climate Informatics.Cham; Springer International Publishing. 2015: 13-21.

[23]NOOTEBOOM P D. Using Network Theory and Machine Learning to predict El Ni?o[J]. 2017,

[24]TROUP A J. The Southern Oscillation Index[J]. Quarterly Journal of the Royal Meteorological Society, 1965,91(390): 490-506.

[25]BAMSTON A G, CHELLIAH M, GOLDENBERG S B. Documentation of a highly ENSO‐related sst region in the equatorial pacific: Research note[J]. Atmosphereocean, 1997, 35(3): 367-383.

[26]CHIODI A M, HARRISON D E. El Ni?o Impacts on Seasonal U.S. Atmospheric Circulation, Temperature, and Precipitation Anomalies: The OLR-Event Perspective*[J].Journal of Climate, 2013, 26(3): 822-837.

[27]LECUN Y, BOSER B E, DENKER J S, HENDERSON D,HOWARD R E, HUBBARD W. Backpropagation Applied to Handwritten Zip Code[J]. Neural Computation, 1989,

[28]LECUN Y. Generalization and Network Design Strategies;proceedings of the Connectionism in Perspective, F, 1989[C].

[29]SHI X, CHEN Z, WANG H, YEUNG D Y, WONG W K,WOO W C. Convolutional LSTM Network: A Machine Learning Approach for Precipitation Nowcasting[J]. 2015,9199: 802-810.

[30]HUBEL D H, WIESEL T N. Receptive fields, binocular interaction and functional architecture in the cat's visual cortex[J]. The Journal of Physiology, 1962, 160(1): 106-154.

[31]FUKUSHIMA K. Neocognitron: A self-organizing neural network model for a mechanism of pattern recognition unaffected by shift in position[J]. Biological Cybernetics,1980, 36(4): 193-202.

[32]LECUN Y, BOTTOU L, BENGIO Y, HAFFNER P. Gradient-based learning applied to document recognition[J]. Proceedings of the IEEE, 1998, 86(11):2278-2324.

[33]KRIZHEVSKY A, SUTSKEVER I, HINTON G E.ImageNet classification with deep convolutional neural networks; proceedings of the International Conference on Neural Information Processing Systems, F, 2012 [C].1097-1105.

[34]KALCHBRENNER N, GREFENSTETTE E, BLUNSOM P. A Convolutional Neural Network for Modelling Sentences[J]. Eprint Arxiv, 2014, 1

[35]CHEN Y, XU L, LIU K, ZENG D, ZHAO J. Event Extraction via Dynamic Multi-Pooling Convolutional Neural Networks; proceedings of the The Meeting of the Association for Computational Linguistics, F, 2015 [C].

[36]KIM Y. Convolutional Neural Networks for Sentence Classi fi cation[J]. Eprint Arxiv, 2014,

[37]REYNOLDS R W, SMITH T M, LIU C, CHELTON D B, CASEY K S, SCHLAX M G. Daily High-Resolution-Blended Analyses for Sea Surface Temperature[J]. Journal of Climate, 2007, 20(22): 5473-5496.

[38]DECKER M, BRUNKE M A, WANG Z, SAKAGUCHI K, ZENG X, BOSILOVICH M G. Evaluation of the Reanalysis Products from GSFC, NCEP, and ECMWF Using Flux Tower Observations[J]. Journal of Climate,2012, 25(6): 1916-1944.

[39]SAHA S, MOORTHI S, PAN H, WU X, WANG J,NADIGA S, TRIPP P, KISTLER R E, WOOLLEN J S,BEHRINGER D. The NCEP Climate Forecast System Reanalysis[J]. Bulletin of the American Meteorological Society, 2010, 91(8): 1015-1057.

[40]SAHA S, MOORTHI S, WU X, WANG J, NADIGA S, TRIPP P, BEHRINGER D W, HOU Y, CHUANG H, IREDELL M. The NCEP Climate Forecast System Version 2[J]. Journal of Climate, 2014, 27(6): 2185-2208.