一種基于時間序列分析的赤潮預測新方法研究
——以浙江海域為例

徐麗麗，余駿，高鑫鑫 ，車助鎂，何雯 ，邱婷

（1. 國家海洋局東海預報中心，上海200081；2. 海洋生態(tài)監(jiān)測與修復技術重點實驗室，上海200081；3. 浙江省海洋監(jiān)測預報中心，浙江杭州211800）

1 引言

赤潮又稱紅潮，國際上也稱為“有害藻華”或“紅色幽靈”。赤潮已成為當今世界普遍關注的海洋生態(tài)問題[1]。我國于2005 年制定了《赤潮海洋災害應急預案》并實施至今，赤潮災害預測及評估工作被列為政府工作報告內容之一；各級海洋減災相關部門和預報業(yè)務機構建立了赤潮年度預測和月及周會商業(yè)務制度，作為我國政府赤潮災害應急決策和處置的工作依據和技術支撐。我國浙江海域受南下長江沖淡水、北上臺灣暖流、錢塘江和甌江等徑流注入以及外海上升流系的綜合影響，其環(huán)境極有利于赤潮生物的生長繁殖。國內最早的赤潮記錄是1933 年費鴻年記載的浙江鎮(zhèn)海至臺州-石浦一帶的夜光藻和骨條藻赤潮[2]。根據1989——2019年《浙江省海洋環(huán)境公報》和相關文獻資料[3]，浙江海域在1981——2019 年間共記錄到赤潮事件744 次，累計面積超過1.74×106km2。尤其進入21 世紀以來，超過上千平方公里的浙江海域內，大型赤潮和有毒赤潮都呈明顯增長趨勢。這不但對人民生命健康與海洋生態(tài)環(huán)境造成了極大威脅，同時也成為浙江省海洋經濟高質量和可持續(xù)發(fā)展的一個重要制約因素，因此引起我國學者和各級政府的高度關注。

至今，學術界已對赤潮發(fā)生過程的模擬預報和赤潮災害趨勢或概率預測等做了大量研究。但因赤潮爆發(fā)原因復雜，各藻種及種間發(fā)展的生態(tài)機制尚不完全清楚，現(xiàn)有的預測研究從開始的主要基于氣象水文等因子的定性分析和條件預報[4-5]，逐漸發(fā)展到利用連續(xù)赤潮現(xiàn)場監(jiān)測或遙感數(shù)據[6-8]、基于主成分和多元回歸等統(tǒng)計模型做預測[9-10]。但在實際業(yè)務應用中，這些方法受限于樣本數(shù)據、變量因子的敏感性等因素，從而影響預報精度乃至難以做出有效預測。隨著監(jiān)測數(shù)據的豐富和某些赤潮種生態(tài)規(guī)律機制研究的深入，部分學者開始利用物理-化學-生物耦合的生態(tài)動力學數(shù)值模型進行赤潮模擬研究[11-17]。如夏綜萬等[16]考慮生物動力學和環(huán)境動力學因素，建立了大鵬灣夜光藻赤潮生態(tài)仿真模型；李雁賓[17]對長江口及鄰近海域季節(jié)性的赤潮生消過程控制機理進行了研究。生態(tài)動力學數(shù)值模型可以對機理明確、變量初始值及邊界條件來源詳盡且連續(xù)的局部海域做某個具體赤潮過程的模擬預測，但在以年為尺度的業(yè)務化預測上的應用效果有限。近年來，以神經網絡和深度學習為代表的大數(shù)據預報技術飛速發(fā)展，其在非線性模式識別方面具有獨特的信息處理和解算能力，非常適用于赤潮這種機制尚不清楚的高維非線性系統(tǒng)[18-19]，但其對數(shù)據要求量大、質量嚴格且建模過程復雜。以上是赤潮預測研究取得的眾多成果，在業(yè)務化應用中也取得一定成效。但由于赤潮爆發(fā)受到水文氣象條件、海水理化因子變化以及船舶帶來的外來浮游物種入侵等眾多因素的影響，加之赤潮生態(tài)系統(tǒng)各因子間表現(xiàn)出的高度非線性和不確定性，以及連續(xù)監(jiān)測數(shù)據獲取困難等問題，目前在一線業(yè)務預報機構中能真正使用，且滿足實際業(yè)務工作需求的赤潮預報工具箱的選擇還不夠多。尤其是針對東海區(qū)浙江海域的業(yè)務化赤潮預測研究主要以定性分析為主[20-23]，業(yè)務一線應用的定量預測方法基本空白。

為滿足自然資源部東海分局及浙江省監(jiān)測預報中心赤潮年度預測和月會商的業(yè)務需求，并規(guī)避數(shù)據限制和機理研究要求，本文擬基于1981——2018年浙江海域赤潮月發(fā)生次數(shù)構建時間序列；考慮赤潮長期的年代變化特征又兼顧其季節(jié)性生態(tài)變化規(guī)律，提出一種基于自回歸移動平均（Auto Regressive Integrated Moving Average，ARIMA）模型的時間序列模型，以及易于實現(xiàn)且便于使用的赤潮預測新方法，為浙江海域的赤潮災害年度預測、業(yè)務化赤潮預警報和災害評估工作提供新的技術工具。

2 數(shù)據與方法

2.1 數(shù)據來源

赤潮原始數(shù)據來源于自然資源部海洋預警監(jiān)測司發(fā)布的1989——2019 年《中國海洋災害公報》及浙江省海洋監(jiān)測預報中心提供的1981——2019 年的監(jiān)測數(shù)據。具體參數(shù)包含赤潮發(fā)生時間、發(fā)生海域、經緯度、最大記錄面積、分布形態(tài)、優(yōu)勢藻種、密度和水色，時間跨度為1981——2019年。因赤潮發(fā)生具有年際、季節(jié)和月等不同時間尺度變化，如果直接利用原始時間序列數(shù)據構造模型，容易因非平穩(wěn)特性產生虛假回歸，因此整理并建立38a 浙江赤潮發(fā)生頻率的時間序列后，先對其進行平穩(wěn)化處理。

2.2 分析方法

時間序列分析法是根據一組相依有序的離散數(shù)據，建立反映時間序列中所包含的動態(tài)依存關系的數(shù)學模型，并進行未來狀態(tài)預測[24]。對于非平穩(wěn)時間序列，主要運用ARIMA 模型，亦稱Box-Jenkins模型。指定3 個參數(shù)，即描述自回歸階數(shù)（p）、差分次數(shù)（d）和移動平均階數(shù)（q），模型通常被寫作ARIMA（p，d，q）。其數(shù)學表達式為：

式中：Zt為原序列；at為白噪聲序列，是一列相互之間無關、均值為0、誤差的方差為σ2的隨機變量序列；B 為后移算子即BZt= Zt-1；φp為自回歸算子，φp( B ) = (1 - φ1B - …- φPBP)，P 為模型的自回歸除數(shù)；θq為移動平均算子，θq( )B = (1 -θ1B - …- θqBq)，q 為模型的移動平均階數(shù)；θ0為參數(shù)，θ0= μ(1 - φ1- φ2- …- φp)，μ為平均數(shù)。

建立模型前需根據時間序列的特性（平穩(wěn)性、非平穩(wěn)性和季節(jié)性）確定建模類型。若序列非平穩(wěn)且有季節(jié)性，則模型函數(shù)被記作ARIMA（p，d，q）（P，D，Q）s。它可以用于分析不僅含有季節(jié)性成分、還混有非季節(jié)性成分的時間序列資料。其中（p，d，q）和（P，D，Q）分別為非季節(jié)性和季節(jié)性自回歸（Auto Regressive，AR）、差 分（I）和 移 動平 均（Moving Average，MA）的階數(shù)，s代表季節(jié)周期。本文中模型的原始序列平穩(wěn)化、模型參數(shù)估計、模型診斷和預測均借助SPSS 25.0 統(tǒng)計分析軟件，采用編程法（語句）分析處理。建模過程見圖1。分3 個關鍵階段[25]：（1）模型參數(shù)的確定。利用自相關和偏自相關分析時間序列原始數(shù)據的隨機性、平穩(wěn)性和季節(jié)性，初步確定模型參數(shù)p、d、q 及P、D、Q 的取值；（2）模型參數(shù)的檢驗。首先采用Box-Ljung 檢驗模型白噪聲判斷模型的擬合優(yōu)度；若有兩個或兩個以上的模型通過檢驗，則根據貝葉斯信息準則（Bayesian Information Criterion，BIC）判斷模型類型和最優(yōu)階次[26]，具體算法見式（2）和（3）；（3）預測應用。通過對比模型預測值與實際值的差值，評價模型預測的準確性。

式中：N為樣本個數(shù)，k為參數(shù)估計的數(shù)量，σ2為誤差的方差。在ARIMA 模型中取BIC 最小值作為評價指標確定p、q以及P、Q值。

圖1 時間序列模型建模具體流程

3 分析與討論

3.1 序列平穩(wěn)化處理

根據1981——2019 年共744 起赤潮發(fā)生事件可知（見圖2），浙江海域赤潮發(fā)生次數(shù)具有明顯的“低頻發(fā)生-爆發(fā)增長-緩降波動”3 段式年際變化特征。20世紀為低頻發(fā)生階段：其中80年代發(fā)生頻次寥寥無幾，年均兩次；90 年代尤其是前期略有增加，年均4.5 次。進入21 世紀后，赤潮發(fā)生次數(shù)呈現(xiàn)爆發(fā)性增長，高達45 次/a；2003 年到達了歷史峰值79 次后逐漸回落。2010 年后進入第三階段，發(fā)生次數(shù)緩降后呈穩(wěn)定波動狀態(tài)，比前10 a明顯減小，年發(fā)生數(shù)維持在18次/a。

從季節(jié)分布來看，一年四季皆有赤潮發(fā)生，但集中爆發(fā)于春夏兩季（發(fā)生次數(shù)占97.5%）。由圖3可知，除11、12月以外，其余各月份均有赤潮發(fā)生且主要在4——8 月。其中5 月發(fā)生最多（287 次），其次為6月（188次），分別占全年的46.7%和29.7%；發(fā)生最少月份為10 月，共發(fā)生4 次，約占0.3%。從持續(xù)天數(shù)來看，浙江海域持續(xù)天數(shù)1～3 d 的短期赤潮出現(xiàn)次數(shù)最多，占64%；20 d以上的超長周期赤潮事件共出現(xiàn)11 次，占1.8%，其中最長持續(xù)天數(shù)為31 d。值得注意的是兩次超長持續(xù)時間的赤潮事件發(fā)生時 間 段分別為2011 年2 月9 日——3 月7 日 和2017 年2 月7 月——3 月9 日，均發(fā)生在冬季的象山港港底海域。相關研究已表明[27]，象山港在電廠建成前冬季平均水溫約為8～9 ℃。2005 年底寧海國華電廠和2006 年底大唐烏沙電廠相繼投產后，溫排水使得附近海域的水溫升高，4 ℃溫升包絡線范圍不斷擴大使水溫達到15 ℃左右的可能性大大增加，而且冬春季溫排水的熱效應使浮游植物量增加約5%。因此推測由于溫排水的熱效應導致生物量增加，從而使得象山港海域冬末春初的赤潮呈現(xiàn)低溫期爆發(fā)、持續(xù)時間長、影響面積小等特點。因此，必須重視特殊海域“低溫期”的赤潮事件，以防在赤潮預警報業(yè)務中“漏報事件”發(fā)生。

圖2 近40 a浙江海域赤潮逐月發(fā)生次數(shù)時序圖

由圖2、3 可見，浙江海域赤潮發(fā)生次數(shù)既存在明顯的年際變化，同時也存在典型的季節(jié)波動，因此導致基于赤潮發(fā)生次數(shù)逐月原始數(shù)據建立的時間序列非平穩(wěn)且具有季節(jié)性。當變量不平穩(wěn)時，若直接構造ARIMA 模型容易產生虛假回歸，因此須先對原始數(shù)據進行一階普通差分。差分后雖然已沒有上升和下降趨勢，但是隨著時間的增加周期性一直存在，因此還需繼續(xù)做季節(jié)差分。圖4顯示，一階季節(jié)差分后，序列的長期趨勢和季節(jié)性趨勢基本消失，數(shù)值圍繞0 上下隨機波動。自相關圖呈現(xiàn)逐漸衰減的趨勢，自相關系數(shù)能夠趨于0，且延遲16階后在0 值附近波動，統(tǒng)計量的相伴概率小于0.05，通過單位根檢驗（ADF）[28]，故可認為處理后的時間序列平穩(wěn)，符合建模的條件。

圖3 近40 a浙江海域赤潮月發(fā)生次數(shù)分布特征圖

圖4 近40 a浙江海域逐月赤潮發(fā)生數(shù)差分序列圖

3.2 模型擬合及參數(shù)確定

上文處理后的平穩(wěn)時間序列，可依據Box-Jenkins的理論方法直接構建ARIMA模型。

根據赤潮發(fā)生次數(shù)時間序列的差分情況及季節(jié)性特征，可以初步判別該時間序列是周期為12的復合季節(jié)模型ARIMA（p，d，q）（P，D，Q）s，故S 取值為12；差分后序列自相關系數(shù)（Auto -Correlation Function，ACF）和偏相關系數(shù)（PArtial Correlation Function，PACF）[29]表明（見圖5），ACF和PACF 均在P 參數(shù)大于l 后驟減，故初步判斷連續(xù)模型為RIMA（1，l，1）。季節(jié)模型的參數(shù)P、Q判斷較復雜，一般情況下超過二階的情況很少見，可以分別取0、1、2 由低階到高階逐個實驗。本文基于SPSS 25.0 統(tǒng)計軟件“時間序列預測”模塊中的專家建模器，通過極大似然法[30]進行估計，初步擬定模型參數(shù)，剔除系數(shù)不顯著的模型，并對剩下模型的殘差進行Ljung-Box非線性檢驗[31]，選出P＞0.05 的模型；再從選出的模型中遵從BIC 準則[26]篩選出最優(yōu)模型。由表1 結果可知，Ljung-Box 的Q 檢驗顯示殘差目前并未違反白噪聲的假設，也沒有出現(xiàn)離群值，選取的赤潮發(fā)生次數(shù)時間序列的最優(yōu)預測模型為ARIMA（1，1，1）（1，1，0）12，擬合優(yōu)度系數(shù)為0.68，結果在可接受范圍內。圖5 可見模型殘差的ACF 和PACF 均≤0.5，殘差序列各數(shù)值間沒有相關性，這說明建立的預測模型已充分提取了序列信息，是合適且可信的。

3.3 模型檢驗及預測討論

一般情況下，為評估預測模型的穩(wěn)定性和適應性，會選擇擬合優(yōu)度、平均絕對誤差和相對誤差來評價模型的整體擬合度[32-33]。因為赤潮逐月發(fā)生次數(shù)樣本的特殊性，大多數(shù)樣本為0值或數(shù)值較小，從而導致絕對誤差和相對誤差較大，但在實際業(yè)務中卻屬于可接受范圍，因此本文采用絕對誤差值范圍出現(xiàn)概率作為評價指標?；谏衔慕⒌淖顑?yōu)模型ARIMA（1，1，1）（1，1，0）12，對1981 年1 月——2018年12月的赤潮逐月發(fā)生次數(shù)進行模擬計算，模型擬合優(yōu)度系數(shù)為0.68，模擬值與實測值較吻合（見圖6）。

圖5 殘差的ACF和PACF圖

表1 浙江海域赤潮發(fā)生次數(shù)ARIMA預測模型的相關參數(shù)

ARIMA 模型中，數(shù)據樣本的時間尺度以月為單位進行計算。744 次赤潮發(fā)生事件按照實際發(fā)生時間歸入各月中，形成ARIMA 模型的457 個樣本。由表2 可知，在457 個樣本中，絕對誤差控制在兩次以內的樣本共382 個，占總數(shù)的84%；其中絕對誤差為0 次即模擬值與實測值完全吻合的月份有292個，占比64%。但也存在個別月份絕對誤差較大的情況，主要出現(xiàn)在2000——2005年浙江海域赤潮爆發(fā)性增長階段，尤其是2004 年5 月赤潮發(fā)生次數(shù)爆發(fā)性增長至峰值59 次，ARIMA 模型未能準確擬合這種超歷史極值的小概率異常情況。

利用上文建立的模型對浙江海域2019 年赤潮逐月發(fā)生次數(shù)進行預測（見圖7）。2019年浙江海域赤潮實際發(fā)生次數(shù)為22 次，預測值為19 次，赤潮年發(fā)生次數(shù)的相對誤差為14%。其中4 月、5 月和8 月的絕對誤差都控制在兩次以內，其余月份預測值與實測值完全吻合（見表3）。 可見本文建立的ARIMA 模型能夠較準確地進行赤潮發(fā)生次數(shù)的年度趨勢預測。

圖6 赤潮逐月發(fā)生數(shù)ARIMA模型預測值與實測值對比

表2 ARIMA模型模擬值與實測值的絕對誤差統(tǒng)計

圖7 浙江海域2019年赤潮逐月發(fā)生次數(shù)預測

表3 浙江海域2019年赤潮逐月發(fā)生次數(shù)預測值與實測值的絕對誤差統(tǒng)計

4 結論

考慮赤潮系統(tǒng)的高維非線性機制和年度預測業(yè)務需求，本文利用近40 a 浙江海域赤潮逐月發(fā)生次數(shù)的時間序列，分析其在年際、季和月3個不同時間尺度的變化特征；通過對原始序列的差分處理，基于時間序列分析方法建立了ARIMA 預測模型，并對2019 年浙江海域赤潮逐月發(fā)生次數(shù)進行后報和檢驗。結果如下：

（1）浙江海域赤潮發(fā)生次數(shù)存在明顯年際變化特征，經歷了“低頻發(fā)生-爆發(fā)性增長-緩降后波動”3個明顯的階段，其中2000 年和2010 年為轉折點。同時也存在典型的季節(jié)波動，97.5% 集中爆發(fā)在春、夏季，其中每年的5 月發(fā)生最多，共298 次，占全年的46.7% 。從持續(xù)天數(shù)來看，1～3 d 的短期赤潮出現(xiàn)次數(shù)最多，占64%。

（2）通過對原始序列的差分處理和參數(shù)檢驗，最終建立ARIMA（1，1，1）（1，1，0）12模型且殘差通過白噪聲檢驗，擬合優(yōu)度系數(shù)為0.68，絕對誤差控制在兩次以內的樣本占總數(shù)的84%。用此模型對浙江海域2019年赤潮發(fā)生次數(shù)進行預測，后報檢驗顯示年發(fā)生總次數(shù)相對誤差為14%，各月絕對誤差均在兩次以內，預測結果與實際較吻合。

（3）時間序列預測法需要的僅是序列本身的歷史數(shù)據，對機理研究沒有高要求，在赤潮年度預測等業(yè)務中具有簡易實用且經濟性好的優(yōu)勢；同時在具備長時間序列數(shù)據的基礎上，其精度在業(yè)務應用可接受范圍內。在今后的相關研究中，可探索線性與非線性模型的最優(yōu)組合模型，如ARIMA 與非線性自回歸神經網絡模型組合，加強赤潮相關影響因子的收集并納入時間序列模型中。這樣既保留模型本身的線性預測能力，又彌補了其在非線性預測方面的不足，以期提高模型預測精度。

（4）本文所建立的ARIMA 模型具有簡捷經濟實用、業(yè)務性強的特點。但同時值得注意的是，因其僅考慮時間序列上的依存性和隨機波動的干擾性，對于超歷史極值會有無數(shù)據依存導致擬合不佳的情況。如2004 年5 月出現(xiàn)歷史第一極值“59 次”，分析其原因主要是海水養(yǎng)殖面積、沿海城市集聚發(fā)展導致近岸海域4類及劣4類海水比例在2004年達到峰值，同期海溫也處于快速上升期。在營養(yǎng)鹽充沛供給、水文氣象條件適宜等多種因素共同作用下導致了赤潮頻繁爆發(fā)，出現(xiàn)超歷史極值。因此，下一步研究中需將時間序列模型的預測結果與其他關聯(lián)因素相結合開展綜合統(tǒng)計分析，進一步改進超歷史極值情況的預測精度。