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

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

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

1 引言

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

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

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

2 數據與方法

2.1 數據來源

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

2.2 分析方法

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

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

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

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

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

3 分析與討論

3.1 序列平穩化處理

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

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

圖2 近40 a浙江海域赤潮逐月發生次數時序圖

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

圖3 近40 a浙江海域赤潮月發生次數分布特征圖

圖4 近40 a浙江海域逐月赤潮發生數差分序列圖

3.2 模型擬合及參數確定

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

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

3.3 模型檢驗及預測討論

一般情況下，為評估預測模型的穩定性和適應性，會選擇擬合優度、平均絕對誤差和相對誤差來評價模型的整體擬合度[32-33]。因為赤潮逐月發生次數樣本的特殊性，大多數樣本為0值或數值較小，從而導致絕對誤差和相對誤差較大，但在實際業務中卻屬于可接受范圍，因此本文采用絕對誤差值范圍出現概率作為評價指標。基于上文建立的最優模型ARIMA（1，1，1）（1，1，0）12，對1981 年1 月——2018年12月的赤潮逐月發生次數進行模擬計算，模型擬合優度系數為0.68，模擬值與實測值較吻合（見圖6）。

圖5 殘差的ACF和PACF圖

表1 浙江海域赤潮發生次數ARIMA預測模型的相關參數

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

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

圖6 赤潮逐月發生數ARIMA模型預測值與實測值對比

表2 ARIMA模型模擬值與實測值的絕對誤差統計

圖7 浙江海域2019年赤潮逐月發生次數預測

表3 浙江海域2019年赤潮逐月發生次數預測值與實測值的絕對誤差統計

4 結論

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

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

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

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

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