趙 健，蔡瑞陽，孫偉富

（1.中國石油大學（華東）海洋與空間信息學院，山東 青島 266580；2.江蘇滿運物流信息有限公司，江蘇 南京 210012；3.自然資源部第一海洋研究所，山東 青島 266061）

氣候變暖背景下全球平均海平面呈持續(xù)上升趨勢，給人類社會的生存和發(fā)展帶來嚴重挑戰(zhàn)，已成為當今國際社會普遍關注的全球性熱點問題［1］。海平面上升是一個緩慢但持續(xù)的過程，其長期累積的結(jié)果將對沿岸地區(qū)構成嚴重威脅［2］。1970年―2015年，全球平均海平面上升了0.93 m［3］。近40年來中國沿海海平面變化總體呈波動上升趨勢，1980年―2021年中國沿海海平面上升速率為每年3.4 mm，遠高于每年2.0 mm的全球平均上升速率［4］。

我國約70%以上的大中城市和一半以上的人口集中分布在東部沿海地區(qū)，該區(qū)域人口密集、經(jīng)濟發(fā)達，但由于地勢較低和地質(zhì)沙質(zhì)等原因，該區(qū)域也是受海平面上升影響最為嚴重的脆弱區(qū)域［5］。王龍等［6］利用衛(wèi)星高度計數(shù)據(jù)研究東海海平面變化的線性趨勢，趙健等［7］提出一種反向傳播（BP）神經(jīng)網(wǎng)絡組合模型研究東海海域海平面長期變化趨勢，Wang等［8］構建了多模態(tài)融合的神經(jīng)網(wǎng)絡模型以預測南海海域海平面變化。目前，針對我國近海海域海平面變化的研究大多集中于海平面信號變化較為豐富的南海以及東海海域，中高緯度黃渤海海域的相關研究則較少。

海平面變化預測方法主要有數(shù)理統(tǒng)計分析和氣候模式預測。氣候模式預測多用于全球及大范圍預測，對區(qū)域海平面變化的預測能力有限；數(shù)理統(tǒng)計分析則是將海平面變化分為確定性趨勢項、確定性周期項、剩余分量和白噪聲序列等，再逐個對原始序列進行擬合和外推，多受制于歷史觀測數(shù)據(jù)的長度、質(zhì)量等因素［9］。隨著人工智能技術的快速發(fā)展，各種深度學習算法不斷更新，許多研究者嘗試利用神經(jīng)網(wǎng)絡實現(xiàn)長時間序列的分類預測［10-12］。作為深度學習算法的典型代表，長短期記憶（long short-term memory，LSTM）網(wǎng) 絡 的 應 用 也 日 益 成 熟［13］。LSTM網(wǎng)絡是一種時間循環(huán)神經(jīng)網(wǎng)絡，能夠很好地針對序列數(shù)據(jù)長期依賴問題進行建模［14］。作為一種數(shù)字信號處理技術，奇異譜分析（SSA）能從時間序列中提取觀測數(shù)據(jù)的非線性趨勢，特別適合分析和預測有周期震蕩的時間序列數(shù)據(jù)［15-16］。

為進一步加強對中國近海中高緯度海域海平面變化的研究，以便全面掌握中國沿海海平面變化規(guī)律，選取中國黃海海域作為研究區(qū)域，利用1993年―2015年的海平面異常（SLA）時空序列，將SSA與LSTM網(wǎng)絡構建的組合模型（下文簡稱SSA-LSTM組合模型）應用到黃海海域海平面變化多尺度預測中，以驗證SSA-LSTM組合模型的有效性。最后，利用1993年―2015年的海平面異常數(shù)據(jù)，基于SSA-LSTM組合模型預測2016年―2025年黃海海域海平面異常變化趨勢。

1 研究區(qū)域與數(shù)據(jù)

黃海是太平洋西部的一個邊緣海（見圖1），位于中國大陸與朝鮮半島之間（122°E―127°E，37°N―40°N以及119°E―127°E，31N°―37°N）。黃海是以旅順口至蓬萊的連線為北邊界、以濟州島至長江口的連線為南邊界的半封閉陸架淺海，平均水深44 m。黃海由山東半島成山角和朝鮮長山串之間的連線分為北黃海和南黃海兩部分，北黃海東西兩側(cè)被陸地包圍，南黃海東側(cè)通過對馬海峽與日本海相通［17］。受太陽輻射、降雨、入海徑流及風應力等因素的影響，黃海平均海平面振幅變化較大，季節(jié)性差異明顯［18-19］。

圖1 研究區(qū)域Fig.1 Study area

海洋氣候數(shù)據(jù)集（CDRs）是“可用于確定氣候變率和氣候變化，具有足夠時間長度、一致性和連續(xù)性的長時間序列數(shù)據(jù)，是開展全球氣候變化研究的有效數(shù)據(jù)集”［20-21］。自然資源部第一海洋研究所融合我國自主發(fā)射衛(wèi)星的遙感數(shù)據(jù)，形成了國內(nèi)首套具有自主知識產(chǎn)權、長時間序列、高時空分辨率和連續(xù)的全球海洋CDRs，為海洋環(huán)境變化監(jiān)測以及全球氣候變化研究提供了科學依據(jù)［21］。利用CDRs中的海平面異常數(shù)據(jù)進行黃海海域海平面變化預測研究，海平面異常數(shù)據(jù)的時間范圍為1993年―2015年，時間分辨率為1個月，空間分辨率為25 km，經(jīng)與驗潮站數(shù)據(jù)考核比較，精度指標為每年0.2 mm。區(qū)域平均后的黃海海域1993年―2015年海平面異常數(shù)據(jù)時間序列如圖2所示，線性擬合可得該時間段黃海海域海平面上升速率約為每年3.67 mm，與前人研究結(jié)果一致［22］。

圖2 區(qū)域平均后的黃海海域1993年―2015年海平面異常時間序列Fig.2 SLA time series of the Yellow Sea from 1993 to 2015 after regional average

2 原理與方法

2.1 奇異譜分析

奇異譜分析是在K-L變換分解理論基礎上發(fā)展起來的，與經(jīng)驗正交函數(shù)有關［23］。奇異譜分析構造了一維非線性時間序列的軌跡矩陣，可以分解和重建該矩陣以提取原始時間序列的各個組成部分，如長期趨勢、周期項或噪聲［24］。分解和重建是互補的階段，先將具有顯著震蕩的分量提取出來，再選擇其中有意義的分量進行重建，具體步驟如下［25］：

步驟1嵌入。選擇合適的窗口長度L將原始序列（x1，x2，…，xN）進行滯后排列，得到軌跡矩陣X，如下所示：

式中：N為序列長度。一般L＜N/2，通常為序列周期的最小公倍數(shù)。

步驟2分解。對軌跡矩陣的協(xié)方差矩陣S=XXT進行奇異值分解，得到特征值λ1≥λ2≥…≥λL≥0和對應的特征向量U=(U1，U2，…，UL)，Ui稱為時間正交函數(shù)，記作TEOF。令Vi=則矩陣X的奇異值分解可以寫為

步驟3分組。根據(jù)提取信息的頻率等成分判斷，將X分為m個不相交的組I1，I2，…，Im，代表不同的趨勢成分。以第j組為例，Ij=(i1，…，io)，有

則X表示為

步驟4重構。計算軌跡矩陣X在Um上的投影，如下所示：式中：ai，m為Xi所反映的時間正交函數(shù)在原序列時段上的權重，稱為時間主成分，記作TPC。根據(jù)時間正交函數(shù)TEOF和時間主成分TPC進行對角平均，重建得到xi，k，如下所示：

2.2 LSTM網(wǎng)絡

LSTM網(wǎng)絡是一種特殊的遞歸神經(jīng)網(wǎng)絡［26］，在隱藏層中使用存儲單元代替神經(jīng)單元來解決梯度消失或梯度爆炸問題。LSTM網(wǎng)絡是由每個長短期記憶結(jié)構單元串接而成，包括控制門單元和輸出門單元，控制門單元又是忘記門單元和輸入門單元的加和組成。具體步驟如下［14］：

（1）決定從細胞狀態(tài)丟棄什么信息。通過忘記門讀取ht-1和xt，輸出一系列0到1之間的數(shù)ft，計算式如下所示：

（2）決定往細胞狀態(tài)存放什么信息。首先通過輸入門決定更新信息it，然后通過tanh層創(chuàng)建一個新的候選值向量，計算式如下所示：

（3）將舊細胞狀態(tài)Ct-1與ft相乘，再加上，根據(jù)更新狀態(tài)的程度進行變化，計算式如下所示：

（4）決定輸出什么信息。首先，通過輸出門決定輸出信息Ot，然后通過tanh層將此時的細胞狀態(tài)Ct和輸出信息Ot相乘得到確定輸出部分ht，計算式如下所示：

式中：xt、ht分別為t時刻網(wǎng)絡的輸入與輸出；ft、it、Ot、Ct分別為忘記門、輸入門、輸出門和記憶細胞的輸出；Wf、Wi、WO、WC分別為對應的權重矩陣；bf、bi、bO、bC分別為對應的偏置；σ為Sigmoid激活函數(shù)。

2.3 SSA-LSTM組合模型

海平面變化的非平穩(wěn)性嚴重影響著預測模型的效果，SSA-LSTM組合模型對區(qū)域海平面變化的預測主要包括：①海平面異常序列的奇異譜分解，降低原始序列的復雜度；②分解后時間模態(tài)子序列的LSTM網(wǎng)絡預測，基于LSTM網(wǎng)絡完成模型預測；③預測結(jié)果的重構，利用奇異譜分析的特點實現(xiàn)海平面時間序列、時空序列的雙預測。奇異譜分析可把隨時間變化的海洋要素場分解為空間模態(tài)和時間模態(tài)兩部分：空間模態(tài)包含要素場的空間分布信息，該部分是不隨時間變化的；時間模態(tài)則由空間點的線性組合所構成，稱為主分量［27］。針對時間模態(tài)建立LSTM網(wǎng)絡預測模型，可有效改善LSTM網(wǎng)絡預測帶來的滯后問題，彌補數(shù)據(jù)分解存在的模式混疊現(xiàn)象。時空序列預測與時間序列預測的不同在于：軌跡矩陣的創(chuàng)建，前者將海平面網(wǎng)格數(shù)據(jù)視為分解對象，窗口長度由研究區(qū)域經(jīng)度的跨度決定；預測數(shù)據(jù)的重構，分解后前者選取時間模態(tài)進行LSTM網(wǎng)絡預測，預測結(jié)果再與之前存儲的空間模態(tài)相乘即得到時空矩陣（最終預測結(jié)果），無需進行對角平均來重構序列。

SSA-LSTM組合模型通過改變軌跡矩陣的創(chuàng)建方式，實現(xiàn)網(wǎng)格數(shù)據(jù)時空方向上的分解；引入LSTM網(wǎng)絡，應用于海平面變化預測，針對不同需求可以實現(xiàn)時間序列、時空序列的雙預測，同時可以保證較高的預測精度。

3 結(jié)果與分析

3.1 區(qū)域海平面異常序列分析

（1）海平面時間序列分析

海平面變化的影響因素較多，是非線性、非平穩(wěn)性且多尺度變化的時間序列，分解重構方法能有效降低原始序列的復雜性，從而提高預測模型的精度［28］。在對黃海海域海平面異常時間序列進行奇異值分解時，根據(jù)該區(qū)域年周期信號顯著的特點選擇12為窗口長度來創(chuàng)建軌跡矩陣，并利用奇異譜分析對其進行分解；根據(jù)奇異值差分譜進行篩選并按照各自的主周期進行合并分組，最終得到5個時間模態(tài)子序列（REC1―REC5）和1個殘差序列（REC6），如圖3所示。奇異譜分析后的時間模態(tài)子序列頻域穩(wěn)定，由頻譜分析可知，各子序列的主周期幾乎沒有重疊在同一頻率上，REC1―REC5的主周期分別對應11.96個月、30.56個月、5.978個月、3.986個月和2.957個月，表現(xiàn)為明顯的年際變化、2～3年周期變化、半年變化和季節(jié)變化。REC2還存在137.5個月（約11年）、91.67個月、68.75個月等周期，準11年的周期反映了太陽黑子的影響，4～7年的周期通常認為與厄爾尼諾-南方濤動事件有關［29］。海平面不同周期的振蕩在不同時間段交叉，導致其上升、下降振幅的不穩(wěn)定性。

圖3 黃海海域海平面異常時間序列SSA分解結(jié)果Fig.3 SSA decomposition results of SLA time series in the Yellow Sea

（2）海平面時空序列分析

時空序列數(shù)據(jù)反映了區(qū)域內(nèi)多個隨機變量隨時間不斷變化的趨勢，除了時間和空間特征，時空數(shù)據(jù)還包含屬性特征，因而呈現(xiàn)出多維、動態(tài)和時空相關等復雜性質(zhì)［30］。以黃海海域海平面網(wǎng)格數(shù)據(jù)為分解對象，構建36×32大小的軌跡矩陣，利用奇異譜分析進行分解，獲取空間分布形態(tài)以及對應的時間變化特征。結(jié)果表明，第一模態(tài)和第二模態(tài)的累計貢獻率達到86.69%，如圖4所示，這很好地詮釋了1993年―2015年黃海海域海平面異常數(shù)據(jù)的分布類型。

圖4a為第一模態(tài)空間分布場，其方差貢獻率為82.24%，遠高于其他模態(tài)的貢獻率，是黃海海域海平面異常數(shù)據(jù)的主要空間分布形式?？梢钥闯?，特征值均為正值，表明1993年―2015年黃海海域的海平面變化趨勢具有高度一致性，高值中心主要位于渤黃海分界處。

圖4b為第二模態(tài)空間分布場，其方差貢獻率為4.65%，也是黃海海域海平面異常數(shù)據(jù)的典型空間分布形式。以36°N左右為界，向北為負值區(qū)，向南為正值區(qū)，結(jié)合時間系數(shù)判斷春夏季多為負相關，即春夏兩季海平面異常高度由南向北遞增；秋冬季多為正相關，即秋冬兩季海平面異常距平由南向北遞減。

圖4 黃海海域海平面異常數(shù)據(jù)正交函數(shù)空間分布Fig.4 EOF spatial distribution of SLAs in the Yellow Sea

為了獲取黃海海域海平面異常數(shù)據(jù)隨緯度和時間的變化規(guī)律，求得該區(qū)域1993年―2015年共23年的每月及分緯度的海平面異常數(shù)據(jù)均值，構建海平面異常數(shù)據(jù)均值的時間-緯度分布圖［31］，如圖5所示?？梢钥吹?，黃海海域海平面高度極大值出現(xiàn)在秋季9月―10月，極小值則出現(xiàn)在春季3月―4月。秋冬兩季，緯度越高海平面異常距平值越??；春夏兩季，緯度越高則海平面異常距平值越大；在33°N和39°N左右出現(xiàn)海平面異常極值點。

圖5 黃海海域海平面異常數(shù)據(jù)均值的時間-緯度分布Fig.5 Time-latitude distribution of mean SLAs in the Yellow Sea

3.2 區(qū)域海平面變化預測模型構建

（1）海平面時間序列預測

采用LSTM網(wǎng)絡與SSA-LSTM組合模型對黃海海域1993年―2015年的海平面異常時間序列進行預測。經(jīng)奇異譜分析后得到的各時間模態(tài)子序列有顯著的主周期特性，將主周期特性作為輸入并利用LSTM網(wǎng)絡進行預測。將原始海平面異常數(shù)據(jù)分為兩部分：數(shù)據(jù)總量的前80%作為訓練數(shù)據(jù)（1993年1月至2011年6月），數(shù)據(jù)總量的后20%作為檢驗數(shù)據(jù)（2011年7月至2015年12月），此時Dropout對小數(shù)據(jù)防止過擬合有很好的效果。整組實驗中，優(yōu)化算法均采用了自適應估計的Adam（adaptive moment estimation），而學習率與衰落因子根據(jù)各子序列的不同需要在訓練中不斷進行修正，通過經(jīng)驗調(diào)試以平方誤差最小為評價準則來確定LSTM網(wǎng)絡的最優(yōu)參數(shù)。

為了全面評價SSA-LSTM組合模型的預測效果，對研究區(qū)域海平面異常時間序列分別采用LSTM網(wǎng)絡與SSA-LSTM組合模型2種方法進行分析，2種方法的預測結(jié)果與原始觀測數(shù)據(jù)的比對如圖6所示。為避免與原始觀測數(shù)據(jù)比較時產(chǎn)生偏差，首先對奇異譜分析后的時間模態(tài)子序列進行歸一化處理，再將全部時間模態(tài)子序列的預測結(jié)果疊加重構，最終得到SSA-LSTM組合模型的預測結(jié)果。從圖6可以看出，雖然LSTM網(wǎng)絡預測結(jié)果的峰值和谷值與原始觀測數(shù)據(jù)均較為接近，但由于延遲原因?qū)е抡`差相對較大，易出現(xiàn)極值偏差；利用SSA-LSTM組合模型進行預測時，奇異譜分析能很好地對原始序列進行平穩(wěn)化及降噪處理，模擬出更接近原始序列的趨勢，預測結(jié)果優(yōu)于LSTM網(wǎng)絡預測結(jié)果，與原始觀測數(shù)據(jù)最為吻合。

圖6 SSA-LSTM組合模型與LSTM網(wǎng)絡預測結(jié)果對比Fig.6 Comparison of prediction results between SSA-LSTM combined model and LSTM network

為進一步分析上述2種方法的預測效果，采用相關系數(shù)（R）、平均絕對誤差（eMAE）和均方根誤差（eRMSE）等指標對預測結(jié)果進行定量評價，如表1所示。LSTM網(wǎng)絡在模擬原始序列趨勢時表現(xiàn)出良好的優(yōu)越性，相關系數(shù)R較大，但存在極值預測偏差大的問題，eRMSE也較大。SSA-LSTM組合模型能很好地模擬原始觀測數(shù)據(jù)的趨勢變化，同時極值預測能力顯著提高，預測結(jié)果與原始觀測數(shù)據(jù)最為吻合，預測效果最佳。

表1 SSA-LSTM組合模型與LSTM網(wǎng)絡的預測誤差Tab.1 Prediction error of SSA-LSTM combined model and LSTM network

（2）海平面時空序列預測

在對海平面時空序列進行預測建模時，考慮到空間模態(tài)的基本不變性，即短時間內(nèi)空間模態(tài)穩(wěn)定不可變，故只需對時間模態(tài)進行建模［32］。利用奇異譜分析將研究區(qū)域時空序列分解為空間模態(tài)和時間模態(tài)以及奇異值的乘積，再利用LSTM網(wǎng)絡對時間模態(tài)進行預測，預測結(jié)果如圖7所示。最后，將預測結(jié)果與空間模態(tài)及奇異值進行重建得到時空序列各網(wǎng)格點的最終預測結(jié)果。表2為基于SSA-LSTM組合模型的黃海海域海平面異常數(shù)據(jù)時空預測誤差，在為期4年的海平面異常數(shù)據(jù)時空預測中，以年為單位分別對每個時段的時空預測結(jié)果進行誤差統(tǒng)計。可以發(fā)現(xiàn)，隨著預測時間的增加，時空預測精度也逐漸降低，尤其是第三預測年份（2013年6月至2014年6月）預測結(jié)果存在較大偏差。通過海平面異常原始觀測數(shù)據(jù)發(fā)現(xiàn)，2013年11月黃海海域海平面異常數(shù)據(jù)較2012年同期低約200 mm，造成極值點偏差較大，從而導致預測精度較低。海平面異常時空序列預測精度略低于時間序列預測精度，主要是因為在奇異譜分解重構時，包含信息較少的殘差隨機序列被舍棄，而且在對各時間模態(tài)分解后的重復性預測操作容易造成誤差累積。

圖7 時間模態(tài)預測值與真實值對比Fig.7 Comparison of time mode between predicted and true values

表2 基于SSA-LSTM組合模型的黃海海域海平面異常時空預測誤差Tab.2 Spatio-temporal prediction error of SLAs in the Yellow Sea with SSA-LSTM combined model

將SSA-LSTM組合模型與現(xiàn)有方法進行比較，可以更好地理解該模型的預測能力。Yavuzdo?an等［33］于2021年提出了基于Coplula函數(shù)的預報方法，該方法能夠以高精度反映海平面異常數(shù)據(jù)的空間相關性，誤差介于22～41 mm之間；2021年周瑋辰等［34］基于U-Net網(wǎng)絡通過非線性的方法疊加上下采樣過程中不同尺度的空間特征，實現(xiàn)了均方根誤差為39 mm的預報效果。綜上，SSA-LSTM組合模型仍有優(yōu)異表現(xiàn)，可以作為海平面異常時空預測的有力模型。

3.3 未來10年黃海海域海平面變化預測

基于SSA-LSTM組合模型，利用1993年―2015年的黃海海域海平面異常數(shù)據(jù)，在對黃海海域海平面異常時空序列預測的基礎上再進行區(qū)域空間平均，預測黃海海域未來10年的海平面變化趨勢，結(jié)果如圖8所示。經(jīng)線性擬合，2016年―2025年黃海海域海平面變化速率約為每年3.65±0.79 mm，黃海海域海平面將呈持續(xù)上升趨勢。2021年中國海平面公報研究表明，1980年―2021年黃海沿海海平面上升速率為每年3.2 mm，預計未來30年黃海沿海海平面將上升60～160 mm［4］。利用我國自主研發(fā)的海洋CDRs數(shù)據(jù)分析黃海海域海平面變化趨勢，研究結(jié)果與上述公報結(jié)論基本相符，表明SSA-LSTM組合模型具有很好的應用價值，同時也驗證了海洋CDRs的有效性。

圖8 基于SSA-LSTM組合模型的黃海海域2016年―2025年海平面異常預測結(jié)果Fig.8 Prediction results of SLAs in the Yellow Sea with SSA-LSTM combined model from 2016 to 2025

4 結(jié)論

（1）SSA-LSTM組合模型改善了LSTM網(wǎng)絡的滯后現(xiàn)象，極值問題得到緩解。奇異譜分析解決了常見分解算法中出現(xiàn)的模式混疊現(xiàn)象，在考慮空間信息的情況下使預測結(jié)果更接近原始觀測數(shù)據(jù)，組合模型具有很好的應用價值。

（2）基于我國首套自主研發(fā)的海洋CDRs，發(fā)現(xiàn)黃海海域海平面異常具有高度一致性，并且與緯度、時間存在一定相關性。預計2016年―2025年黃海海域海平面變化呈持續(xù)上升趨勢，上升速率約為每年3.65±0.79 mm，與中國海平面公報研究結(jié)果一致，驗證了SSALSTM組合模型和CDRs數(shù)據(jù)的有效性。

作者貢獻聲明：

趙?。赫撐目傮w設計，論文審閱與修改。

蔡瑞陽：數(shù)據(jù)處理，論文初稿撰寫。

孫偉富：海洋氣候數(shù)據(jù)集數(shù)據(jù)處理與分析。