基于多變量LSTM 神經網絡模型的PDO 指數預測研究

于振龍，許東峰, ，姚志雄, ，楊成浩, ，劉松楠

(1.自然資源部第二海洋研究所 衛星海洋環境動力學國家重點實驗室，浙江 杭州 310012；2.浙江省海洋科學研究院，浙江 杭州 310012；；3.自然資源部海洋空間資源管理技術重點實驗室，浙江 杭州 310012)

1 引言

太平洋年代際振蕩（Pacific Decadal Oscillation,PDO）是北太平洋年代際尺度上最主要的海表面溫度（Sea Surface Temperature,SST）異常模態，PDO 的概念最早由Mantua 等[1]于1997 年首次提出，PDO 作為一種太平洋氣候變化強信號，對太平洋及其周邊地區氣候變化有著重要影響。在PDO 暖相位期間，冬季阿留申低壓加強，北太平洋呈現“馬蹄形”的海表面溫度異常，北太平洋中部偏冷，熱帶中東太平洋和北美沿岸偏暖；而在PDO 冷相位期間情況相反，北太平洋中部偏暖，熱帶中東太平洋和北美沿岸偏冷。PDO 指數被定義為太平洋20°N 以北海表面溫度異常場經過經驗正交函數分析（Empirical Orthogonal Function Analysis,EOF）第一模態的時間系數，其正、負分別代表了PDO 的暖、冷位相。PDO 作為主要的年代際氣候變異之一，對太平洋動力環境、漁業收獲[1-3]和我國的降雨[4-5]、氣溫[6-7]、農作物產量[8]以及南海的鹽度變化[9-10]等有重要影響，所以PDO 的預測就顯得尤為重要。

由于PDO 的復雜性以及全球變暖背景下PDO 的可預測性正逐漸降低[11]，預測PDO 的相位變化對于人們來說仍具有很大的挑戰性。前人做過很多關于PDO 形成機理以及PDO 指數重建的研究：Mantua等[1-2]重點分析了PDO 的時間、空間特征，并定義PDO年代際指數是由PDO 指數5 年滑動平均獲得。Yasuda[12]利用樹木年輪和太平洋珊瑚數據對PDO 指數進行了氣候重建。MacDonald 和Case[13]利用加利福尼亞和阿爾伯塔的柔枝松的水文敏感年輪年表重建了公元993-1996 年的PDO 指數時間序列，同時評估了PDO 強度和周期的長期變動特征。趙傳湖等[14]利用東亞地區高分辨率的氣候代用指標，通過雙重檢驗的逐步回歸方法和小波分析等方法，重建了近400年來的北太平洋PDO 指數，并分析了PDO 周期的演變特征。結果表明，北太平洋PDO 具有準20 年和準50～70 年的變化周期。Chen 等[15]通過模型數據比較，探討了千年尺度PDO 變率，并利用北太平洋副熱帶海溫的東西向差異以及東亞地區與PDO 相關的降水異常重建了PDO 指數。Schneider 和Cornuelle[16]指出，PDO 不是單一模式，而是多個物理過程的組合，PDO 的位相轉變源于Ni?o 3.4 的強迫作用以及阿留申低壓和黑潮-親潮延伸體海區變化的累積，并進一步估算了這些影響因子對PDO 的相對貢獻。與Schneider 和Cornuelle[16]的結論一致，Newman 等[17]分析了驅動PDO 的因素，改進了以前的PDO 經驗模型[18]。Huang 和Wang[19]利用海表面高度、海平面壓力和海冰密集度，通過“年際增量法”將變量的年際增量作為預測對象，實現了對PDO指數的后報。Lu 等[20]利用一種氣候網絡分析的新技術來捕獲PDO 相位轉變的前兆。通過計算PDO空間格局冷區和暖區之間的聯系，在相位轉變前幾年，發現冷暖區域之間的負相關逐漸增強。該信號從1890 年到2000 年捕獲了6 個PDO 相位轉變，提前預警時間平均為（6.5±2.3）年。

PDO 是多個物理過程綜合作用的結果。其中，大氣強迫超前于海表溫度變率，是影響北太平洋溫度變率的主要因素。阿留申低壓通過影響北太平洋地區水汽輸送[21]，間接影響北太平洋溫度，是驅動PDO 變率的主要大氣強迫之一；PDO 相位的改變影響風場[22]和溫度場，進而帶來海表面高度的變化[23]；熱含量（Ocean Heat Content，OHC）是表征海洋熱狀態的一個重要參數[24]，海洋主要的熱量存儲于0～700 m 水層[25-26]。大西洋多年代際振蕩（Atlantic Multi-Decadal Oscillation,AMO）可以通過大氣遙相關影響太平洋副熱帶熱量分布，進而影響PDO[27-28]；北極濤動（Arctic Oscillation,AO）在PDO 低頻變化中起著重要的作用，AO 增強會導致阿留申低壓增強，進而通過北太平洋海氣相互作用影響PDO[29]。格陵蘭海冰可以通過影響AO 來影響PDO 變率，海冰密集度可以用來表征海冰的變化[30]。綜上所述，海平面氣壓、海表面高度、熱含量和海冰密集度均可作為PDO 指數的預報要素進行建模。

隨著人工智能技術的發展，機器學習被應用于各個領域的研究[31-33]，深度神經網絡成為了一種新型的研究手段。對于時間序列的預測應用最多的是循環神經網絡（Recurrent Neural Network，RNN）模型，由于RNN 模型結構比較簡單，使用的反向傳播（Back Propagation Through Time，BPTT）算法不能適應長時間序列數據的訓練，容易產生梯度爆炸以及梯度消失的問題，無法實現長期記憶的效果。因此，需要一個含有存儲單元、存儲記憶功能的神經網絡模型。1997 年長短期記憶（Long Short Term Memory，LSTM）神經網絡模型由Hochreiter 和Schmidhuber[34]提出。LSTM 神經網絡模型是RNN 模型的改進，通過增加記憶單元里的控制門限來解決RNN 模型短期記憶的問題，從而可以高效地處理長時間序列信息，實現對復雜函數的建模。目前，國內外利用LSTM 神經網絡模型對ENSO 預報的研究比較多[35-37]，與其他統計預報模型相比，LSTM 神經網絡模型在長時間序列上具有較強的預報能力。由于年代際變化信號的預測比較困難，部分學者主要利用氣候模型實現對PDO 指數的后報試驗[38-39]，利用深度學習的研究還比較少。本文旨在前人研究基礎上，將海平面氣壓、海表面高度、海冰密集度作為PDO 指數的預報要素，新增0～700 m 水層的熱含量作為預報因子，建立LSTM 神經網絡模型，實現對PDO 指數的預測，為PDO 指數的預測研究提供一種新的手段。

2 研究方法

2.1 LSTM 神經網絡模型

LSTM 神經網絡模型在本質上是RNN 模型的改進版本，其與RNN 模型最大的區別是增加了3 個邏輯控制單元（即控制門，包括輸入門、輸出門和遺忘門）來解決RNN 模型短期記憶問題，通過改變神經網絡的記憶單元與其他部分連接的權值大小來控制信息流的輸入、輸出以及記憶單元的狀態。其中輸入門（input gate，記為it）用來控制信息是否能夠被記錄到記憶單元中；輸出門（output gate,記為ot）控制當前時間節點記憶單元中的信息是否能夠流入到當前隱藏層中；遺忘門（forget gate,記為ft）選擇性遺忘上一時刻記憶單元的信息并累加到當前時刻記憶單元中。LSTM 神經網絡模型結構示意圖如圖1 所示，其中，Xt為當前狀態下數據的輸入，ht為當前隱藏層，σ 和 tanh為當前記憶單元中兩個不同的激活函數，“×”和“＋”分別代表矩陣相乘和矩陣相加運算。

圖1 LSTM 神經網絡模型結構示意圖（來源于文獻[34]）Fig.1 Schematic diagram of LSTM neural network model structure (from reference [34])

在LSTM 神經網絡模型的訓練過程中，首先模型把t時刻的數據特征信息輸入到輸入層中，經過激活函數得到輸出結果。然后將輸出結果、t-1 時刻的記憶單元存儲信息和t-1 時刻隱藏層輸出結果輸入到當前LSTM 神經網絡模型結構處理單元中，通過3 個控制門限的處理后，輸出結果將進入下一輸出層和隱藏層，最后計算反向傳播誤差，更新各個遞歸連接權重的大小。在t時刻LSTM 神經網絡模型每個記憶單元中的數據處理過程如式（1）～（6）所示：

式中，ft、it、ot為3 個控制門限；Wf、Wi、Wo分別代表了其對應的控制門限的遞歸連接權重；bf、bi、bo是通過對應的控制門限誤差值；ht代表t時刻隱藏層輸出值；Ct代表當前神經元記憶狀態的記憶單元，具有隨著訓練不斷存儲、讀取、重置和更新長時間信息的能力。控制門結構由激活函數與點乘組成，公式中的激活函數 σ在0～1 之間取值，tanh 在-1～1 之間取值。點乘決定了記憶單元可以傳送過去的信息量。激活函數的計算公式如式（7）和式（8）所示：

2.2 精度評估

利用PDO 預測值與觀測值時間序列的均方根誤差（RMSE）、平均絕對誤差（MAE）和平均絕對百分比誤差（MAPE）評價模型預測值與觀測值的精度大小，計算公式為

3 數據來源及處理

本文采用多個預報要素建立多變量LSTM 神經網絡模型預測PDO 指數，實驗流程如圖2 所示。為了能夠預測更長時間序列的PDO 指數，選取了1921-2020 年100 年的預報要素數據資料。除了海冰密集度數據以外，其他預報要素都選取了相同的研究區域（20°～60°N，120°E～100°W）。

圖2 實驗流程圖Fig.2 Experimental flowchart

首先，PDO 指數本文采用的是美國國家海洋和大氣管理局（National Oceanic and Atmospheric Administration,NOAA）基于第五代版本ERSST 擴展重建的NCEI PDO 指數數據（https://www.ncdc.noaa.gov/teleconnections/pdo/）；海平面氣壓（Sea Level Pressure,SLP）數據源自哥白尼氣候變化服務（Copernicus Climate Change Service，C3S，https://cds.climate.copernicus.eu/）的CMIP5 氣候模式月平均數據，該數據由NOAA 的GFDL-CM2p1 模型計算而來，水平分辨率為2.5°×1.5°；海表面高度（Sea Surface Height,SSH）數據同樣源自C3S 的CMIP5 模式數據，水平分辨率為1°×1°；海冰密集度（Sea Ice Concentration,SIC）數據源自英國氣象局哈德萊中心的全球海冰和海表面溫度數據集HadISST，水平分辨率為1°×1°。

熱含量長時間序列數據獲取相對困難，本文利用熱含量計算公式獲取1921-2020 年的熱含量時間序列，所采用的數據為哈德萊中心的EN4 再分析溫鹽數據集。該數據集融合了WOD9、GTSPP、ASBO、GDACs 等數據集，并進行了質量控制。水平分辨率為1°×1°，時間分辨率為1 個月，深度共42 層，最深達5 000 m 以上。熱含量計算方法為

式中，OHC 為熱含量，單位：J；S為數據網格的面積大小；ρ為海水的密度；Cp為定壓比熱容；T為海水溫度；h為熱含量計算水深（本文中h取為700 m）。

其次，采用相關性分析計算預報要素研究區域內的所有網格點時間序列與PDO 指數的相關系數，根據相關系數選取4 個預報區域，以下稱“預報因子”，其中，SLP 選取的是20°～30°N，160°E～180°海域范圍；SSH選取45°～55°N，175°E～175°W海域范圍；0～700m的OHC選取25°～40°N，150°～170°W海域范圍；SIC選取75°～85°N，10°W～10°E海域范圍。100年的數據時間尺度會大大降低預報要素與PDO 指數的相關性，預報要素與PDO 指數的相關系數空間分布及預報因子區域范圍如圖3 所示，預報因子在區域內每個數據格點的時間變化與PDO 指數均滿足99%置信水平的相關性。

圖3 1921-2020 年預報要素與PDO 指數的相關系數Fig.3 Correlation coefficients of forecast elements and PDO index from 1921 to 2020

最后，計算預報因子的1921-2020 年的時間序列（異常），由于PDO 是年代際尺度上氣候內部變率，周期通常為20～30 年，本文將PDO 指數數據利用Butterworth 低通濾波器進行了5 年的低通濾波處理，去掉高頻變化信號得到PDO 指數，4 個預報因子的數據處理與PDO 指數相同，這樣不僅能夠讓LSTM 神經網絡模型預測更長時間尺度的PDO 指數序列，同時可以提升模型訓練時的訓練效率。圖4 中從上到下依次是PDO 指數、海平面氣壓異常、海平面高度異常、0～700 m 熱含量異常以及海冰密集度異常數據。

圖4 太平洋年代際振蕩（PDO）指數以及預報因子時間序列Fig.4 Pacific decadal oscillation index and predictor time series

4 模型訓練與調參

本文全部數據采用歸一化處理以消除量綱和數量級帶來的影響，有利于提高預測模型的收斂速度和精度。多變量LSTM 神經網絡模型預測時，預測的精度對模型參數的要求較大,不同的參數組合類型對模型的預測精度影響較為明顯。LSTM 神經網絡模型的性能取決于模型中超參數的設置，超參數是在模型訓練前手動設定的訓練變量，包括激活函數、優化器、批大小、訓練輪數和隱藏層節點數。為了提高模型的訓練效率，本文提前確定不重要的超參數：模型設定訓練輪數為300 輪；學習率設置為0.005；使用自適應矩估計（Adaptive Moment Estimation,Adam）優化器不斷更新訓練網絡權重以使損失最小。由于數據量大且為復雜矩陣，故調用GPU 進行訓練加速。為了防止LSTM 神經網絡模型出現過擬合的情況，通過在全連接層后添加隨機失活（dropdout）層對在正相傳遞或權值更新過程中的LSTM 神經網絡模型記憶單元進行優化，隨機失活率的大小設置為0.2。

關于模型最優參數的確定，本文對隱藏層節點數、批大小以及學習率下降周期3 個超參數進行了組合調參，確保訓練精度的同時也降低了模型訓練的復雜程度。采用控制變量的方法，將隱藏層節點數分別取120層、240 層和360 層，批大小取60、120 和180，學習率下降周期取75、150 和225 輪，共27 種參數組合。每種參數組合在模型建立時都需要通過交叉驗證（Cross Validation，CV），交叉驗證是一種流行的具備魯棒性的模型性能評估技術。本文建立的是關于PDO 指數時間序列的多變量多步LSTM 神經網絡模型，訓練集和測試集在邏輯上存在時間先后順序，采用“遞增時間窗交叉驗證法”對模型進行交叉驗證，方法如圖5 所示。

圖5 遞增時間窗交叉驗證法Fig.5 Incremental time window cross-validation

該交叉驗證方法將PDO 指數數據等分成10 份，考慮到PDO 周期在20 年以上，本文將最小訓練集設置為20 年。首先將前20 年（1921-1940 年）作為模型的訓練集，后10 年（1941-1950 年）作為模型的測試集；將前30 年（1921-1950 年）作為模型的訓練集，后10 年（1951-1960 年）作為模型的測試集；以此類推，直到將前90 年（1921-2010 年）作為模型的訓練集，后10 年（2011-2020 年）作為模型的測試集。然后分別計算測試集中預測值與實際值的相關系數及均方根誤差，模型訓練5 次后取平均結果。最終將交叉驗證的結果取平均作為該參數組合的交叉驗證結果，通過對比不同參數組合的模型平均相關系數以及均方根誤差大小，確定模型的最優參數并利用最優參數建立LSTM 神經網絡模型，模型訓練10 次取平均作為最終預報結果，最終將1921-2020 年的預報因子數據作為輸入，實現對未來10 年的PDO 指數預測。

5 結果分析

多變量LSTM 神經網絡模型調參的過程實際上是模型尋找靶心的過程，經過交叉驗證，根據測試集與訓練集相關系數（R）以及均方根誤差（RMSE）確定模型的最優參數，可以最大程度地保證模型的靶心處于最佳位置。模型調參結果如表1 所示，擬合效果比較優異，所有參數組合測試集平均相關系數都在0.6～0.7 之間。其中，隱藏層節點數為360，批大小為180，學習率下降周期為75 輪的模型的平均相關系數為0.697 2，均方根誤差為0.622 1，選取該參數組合作為模型訓練的最優參數，已在表1 中粗體顯示。

表1 交叉驗證平均結果Table 1 Cross-validation average result

圖6 為根據最優參數建立的LSTM 模型10 次訓練的平均擬合及預測結果。1921-2010 年為模型訓練集擬合結果，2011-2020 年為模型測試集的擬合結果。黑色曲線代表NOAA 提供的PDO 指數觀測值；紅色曲線代表多變量LSTM 神經網絡模型訓練結果（預測值）。PDO 指數預測值與觀測值的變化趨勢基本一致，在峰值和谷值處的差距較大。2013-2019 年的預測值比觀測值偏小，其余時間的預測值比觀測值偏大，總體來看，LSTM 神經網絡模型的擬合效果比較理想；藍色曲線代表LSTM 神經網絡模型對未來10 年（2021-2030 年）的低頻PDO 指數預測結果。預測結果顯示，PDO 未來10 年將一直處于冷位相，北太平洋中部出現SST 正異常，熱帶中東太平洋和北美沿岸出現SST 負異常，PDO 指數出現兩次波動，于2025 年出現最小值。

圖6 LSTM 網絡模型預測值及觀測值Fig.6 LSTM neural network model prediction results and true values

如圖7 所示，本文對比了不同模型對2011-2020 年PDO 指數的預測結果，包括單一變量LSTM 神經網絡模型，統計回歸模型以及多變量支持向量回歸（Support Vector Regression，SVR）模型。其中，單變量LSTM 神經網絡模型是由PDO 指數原始時間序列（來自NOAA）構建，經過與多變量LSTM 神經網絡模型相同的交叉驗證方法確定模型最優參數組合，并訓練10 次取平均作為最終結果；統計回歸模型采用多元回歸，與多變量LSTM 神經網絡模型采用相同的預報因子作為變量；多變量SVR 模型變量的選取與本文采用的多變量LSTM 神經網絡模型相同。結果顯示，多變量LSTM 神經網絡模型擬合效果優于其他模型，特別是在峰值和谷值處更接近實際值。圖8 為不同模型的誤差精度評估結果，多變量LSTM 神經網絡模型RMSE=0.358 0，MAE=0.298 4，MAPE=0.406 9；用PDO 指數原始時間序列構建的單變量LSTM 神經網絡模型RMSE=0.447 2，MAE=0.342 4，MAPE=0.529 1；統計回歸模型RMSE=0.466 1，MAE=0.409 9，MAPE=0.590 1；多變量SVR 模型RMSE=0.484 1，MAE=0.424 2，MAPE=0.629 9。從對比結果可以看出多變量LSTM 神經網絡模型預測效果最好，誤差最小，與實際值擬合效果最好，其次是單變量LSTM 神經網絡模型、統計回歸模型、多變量SVR 模型。

圖7 2011-2020 年不同模型的預測結果Fig.7 Forecast results of different models from 2011 to 2020

圖8 不同模型預測的精度評估結果Fig.8 Accuracy evaluation results of different model predictions

6 結論

利用1921-2020 年的海平面氣壓、海平面高度、熱含量數據以及海冰密集度作為PDO 指數的預報要素，通過空間相關性分析選取PDO 指數預報因子，進一步將預報因子的時間序列數據進行5 年低通濾波處理，通過交叉驗證選取最優參數組合，建立關于PDO 指數預測的多變量LSTM 神經網絡模型，對比分析不同模型2011-2021 年的PDO 指數預測結果，進而利用模型實現了對2021-2030 年的PDO 指數的預測。通過此次研究，本文得到如下結論。

（1）基于多個預報因子數據建立的多變量LSTM神經網絡模型可以有效地實現對PDO 指數的預測，模型經過交叉驗證后平均相關系數和均方根誤差分別為0.70 和0.62，預測結果較為理想。

（2）多變量LSTM 神經網絡模型預測效果要優于多變量SVR 模型、統計回歸模型以及單變量LSTM神經網絡模型，其PDO 指數預測結果誤差小、擬合效果好，可以作為一種新型的預測PDO 指數的手段。

（3）本文預測了未來10 年的低頻PDO 指數結果，結果顯示，PDO 未來10 年將一直處于冷位相，北太平洋中部將出現SST 正異常，熱帶中東太平洋和北美沿岸出現SST 負異常，PDO 指數出現兩次波動，于2025年出現最小值。

PDO 作為年代際氣候變化信號，本文建立的多變量的LSTM 人工神經網絡模型實現了PDO 指數的位相預測。前人的研究[19]雖然能夠有效向后預報1975-2009 年的PDO 指數（相關系數為0.81），但是其算法并不能實現對未來年代際時間尺度上PDO 指數的預報。本研究優勢在于采用機器學習方法，采用更大空間范圍的預報因子，一定程度上克服了預測結果對PDO 的依賴，實現了對未來10 年PDO 指數的預報。然而本研究的交叉驗證試驗結果與實測值相關系數在0.6～0.7 之間，低于Huang 和Wang[19]預測的相關性，產生這種結果的原因可能是其預測對象為PDO 冬季（12 月、1 月、2 月）數據，而本研究選取所有月份的數據；其次，超參數的選取具有一定主觀性，在今后的研究中可以適當增加參數組合的數量來減小這種誤差；最后，長時間尺度理想的預測結果要求模型用更長時間的訓練集，在更長時間尺度上的預測精度有待提高。