劉旭，梁穎祺，王兆毅，李志杰，王崢，季軒梁，何恩業(yè)

（1.國家海洋環(huán)境預報中心，北京 100081；2.北京林業(yè)大學經(jīng)濟管理學院，北京 100083）

1 引言

自2007 年首次在青島監(jiān)測到滸苔規(guī)模性聚集形成綠潮以來[1-2]，滸苔災害已經(jīng)連續(xù)10余年在黃海周期性爆發(fā)。大量滸苔堆積在近岸海域，給海洋生態(tài)環(huán)境造成嚴重威脅[3-4]。研究表明海表面溫度（Sea Surface Temperature，SST）是影響滸苔生長和死亡的主要生態(tài)因子[5-6]。通過實驗室設置溫度梯度可以分析滸苔對溫度的敏感性，吳洪喜等[7]認為滸苔適宜生長的海水溫度為10～30 ℃，在缺乏陽光的干露狀態(tài)下可承受的最高溫度為38 ℃。Taylor等[8]認為滸苔適宜生長的溫度范圍為10～25 ℃，最高生長率溫度范圍為15～20℃，最適萌發(fā)溫度為20 ℃。韓紅賓等[5]認為滸苔孢子/配子萌發(fā)適宜溫度為15～25 ℃，最適萌發(fā)溫度為20 ℃；滸苔孢子/配子放散的溫度為20～30 ℃，最適宜放散溫度為25 ℃。忻丁豪等[9]認為20～26 ℃是滸苔生長發(fā)育的最佳水溫，此溫度區(qū)間下當鹽度為26～32 時生長率可達到36%以上[10]。陳月紅[11]發(fā)現(xiàn)滸苔在25℃組和30℃組中，相對增長率下降趨勢明顯，說明滸苔不具有耐高溫性[5]，高SST極大地抑制了滸苔的生長。

上述研究均表明SST 是滸苔生長的重要限制因子，但不同學者通過對照培養(yǎng)實驗獲取的滸苔最適生長溫度范圍具有明顯差異性。黃海滸苔爆發(fā)是自然海域多種環(huán)境因子共同影響的過程[12]，因此實驗室模擬滸苔生長環(huán)境存在一定局限性[13]。在實證研究方面，李弘毅[13]基于時間序列的GOCI（Geostationary Ocean Color Imager）影像定性分析了2011—2017年滸苔爆發(fā)的環(huán)境驅(qū)動因素，得出5—7 月SST 都處于滸苔最適宜生長的溫度范圍內(nèi)。張廣宗等[14]基于環(huán)境衛(wèi)星HJ-1A/1B 的CCD（Change-Coupled Device）影像和MODIS（MODerate Resolution Imaging Spectroradiometer）影像對2011—2017年南黃海海域滸苔信息進行了提取，結(jié)合ESRL（Earth System Research Laboratroy）海溫數(shù)據(jù)對滸苔生長進行了定性分析，結(jié)果表明每年5—8 月SST逐月升高為滸苔生消過程提供了環(huán)境條件。辛蕾等[15]和劉旭等[16]基于多源衛(wèi)星數(shù)據(jù)將滸苔生命周期分為發(fā)生—發(fā)展—爆發(fā)—衰落—消亡5 個階段，并定性分析了綠潮覆蓋面積隨表層水溫的變化呈現(xiàn)出先上升后下降的趨勢，指出約在22 ℃時覆蓋面積達到最大值。何恩業(yè)等[17]在構(gòu)建黃海滸苔漂移輸運模型時耦合了生長消亡過程的生態(tài)模塊，將海溫因素耦合到黃海滸苔的生長率公式中，設置最適宜溫度區(qū)間為15～20 ℃，生長率為1.15～1.80，臨界死亡溫度為25 ℃，最大死亡率為1.10。對于綠潮和海溫的關(guān)系，盡管學者在室內(nèi)培養(yǎng)實驗和生態(tài)動力學模型等方面開展了定量研究，但目前基于遙感影像資料的實證研究仍處在定性分析階段，定量化研究鮮有報道。為給綠潮災害規(guī)模預測提供借鑒，本文基于2010—2019 年多源衛(wèi)星遙感黃海綠潮資料，分析了綠潮覆蓋面積與SST 間的協(xié)整關(guān)系和Granger 因果關(guān)系，使用協(xié)整模型和傳遞函數(shù)模型對綠潮災害爆發(fā)條件下的規(guī)模變化過程進行模擬和預測。2018—2019 年綠潮覆蓋面積的預測檢驗可驗證該方法的適用性和可靠性。在滸苔的生消機制尚不完全清晰的背景下，基于數(shù)據(jù)驅(qū)動的時間序列法可有效提取歷史信息，為綠潮預測提供技術(shù)支持，為深入挖掘綠潮的生態(tài)影響因素提供參考。

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來源

黃海綠潮覆蓋面積數(shù)據(jù)資料來源于2010—2019年國家衛(wèi)星海洋應用中心的每日業(yè)務化綠潮遙感影像解譯結(jié)果，其中2010—2018年的反演結(jié)果主要基于MODIS 和HY-1B 數(shù)據(jù)，2019 年反演結(jié)果主要基于MODIS和HY-1C數(shù)據(jù)，MODIS和HY-1B分辨率為250 m，HY-1C 分辨率為50 m。表層海溫數(shù)據(jù)來源于日本氣象廳（Japan Meteorological Agency，JMA）發(fā)布的近實時的日變化表層海溫數(shù)據(jù)（Merged satellite and in-situ data Global Daily Sea Surface Temperature，MGDSST，網(wǎng)址：http://ds.data.jma.go.jp/gmd/goos/data/pub/JMA-product/）。該數(shù)據(jù)融合了AVHRR（Advanced Very High Resolution Radiometer）衛(wèi)星紅外傳感器數(shù)據(jù)、WindSat、AWSRE和AWSR-2 微波傳感器數(shù)據(jù)以及來自浮標和船測現(xiàn)場的全球SST 產(chǎn)品，水平分辨率為0.25°，時間分辨率為1 d，其資料已被廣泛應用于海洋科學研究[18]。

由于綠潮遙感受云層遮擋影響較大，難以獲得成像條件較好的日數(shù)據(jù)，因此，按照每年綠潮生消的起止規(guī)律，將研究周期設定為每年的5 月8 日—8月7 日[19]，挑選出成像條件較好的反演結(jié)果，按照每月1—7 日，8—14 日，15—21 日，22 日—月底進行周平均，形成綠潮遙感覆蓋面積周平均時間序列{yt}。根據(jù)當日滸苔覆蓋區(qū)域的經(jīng)緯度范圍解析對應區(qū)域內(nèi)的日均海溫數(shù)據(jù)，按照綠潮覆蓋面積周平均序列的處理方法，合成相應的SST周平均時間序列{xt}。

將2010—2017 年數(shù)據(jù)作為實驗集用于模型模擬，2018—2019 年數(shù)據(jù)作為測試集用于模型預測檢驗，以此評估模型預測精度。

2.2 研究方法

2.2.1 Granger 因果檢驗

采用Granger 因果檢驗，從統(tǒng)計角度識別自變量海溫序列{xt}是否為因變量綠潮遙感覆蓋面積序列{yt}的原因，即SST 變化是否對綠潮覆蓋面積變化有顯著影響。假設{xt}和{yt}是寬平穩(wěn)序列，至少存在一個h= 1,2…使得式（1）成立，說明{xt}序列歷史信息的加入能提高{yt}的預測精度，則稱xt是yt的原因[20]。

式中，It表示在t時刻及以前的所有信息集合。設yt(h|It)為信息It可獲得的h步最小均方誤差預測，則相應的均方預測誤差記為∑y(h|It)，It{xs│s≤t}為除開變量xs在t時刻及以前的信息。

原假設{xt}不是{yt}的Granger 原因，備擇假設{xt}是{yt}的Granger原因。

在備擇假設成立的條件下，構(gòu)造F統(tǒng)計量：

式中，SSEr為有約束條件下的隨機波動，SSEu為無約束條件下的殘差平方和；q為{xt}序列的歷史延遲階數(shù)；n為序列長度；p為{yt}的自回歸階數(shù)。當F統(tǒng)計量大于F1-α(q,n-q-p- 1)時，拒絕原假設，認為{xt}是{yt}的Granger原因。

2.2.2 協(xié)整模型

協(xié)整模型可以有效地衡量序列之間是否具有長期相關(guān)性。假定海溫時間序列為自變量{xt}，綠潮覆蓋面積序列為因變量{yt}，構(gòu)造回歸模型：

如果回歸殘差序列{ut}平穩(wěn)，即?t～I(0)，則相應變量序列{yt}與自變量序列{xt}之間具有協(xié)整關(guān)系。檢驗雙變量協(xié)整關(guān)系時，需滿足xt～I(d)且yt～I(d)。在建模時首先需要檢驗數(shù)據(jù)的平穩(wěn)性以確定序列是否符合建模要求。當變量均為d階單整時，采用EG 檢驗法對回歸殘差進行ADF（Augmented Dickey Fuller）平穩(wěn)性檢驗，假設條件為：

由于{?t}序列是回歸模型式（4）中最小二乘法估計的殘差，ADF 檢驗法的統(tǒng)計量極限分布與觀測樣本有所不同，因此需采用專門用于協(xié)整關(guān)系的ADF-t檢驗臨界值表判斷{?t}是否存在單位根及協(xié)整關(guān)系[20]。

當式（4）中回歸殘差序列{ut}不平穩(wěn)時，說明蘊含著歷史信息之間的相關(guān)性，可以進一步構(gòu)建{?t}的自回歸滑動平均模型（Auto-Regressive and Moving Average model，ARMA）。協(xié)整模型可表示為：

式中，L(B)為{?t}q階移動平均系數(shù)多項式；A(B)為{?t}p階自回歸系數(shù)多項式；?t為白噪聲序列，?t～N(0,σ2)。

2.2.3 傳遞函數(shù)

建立海溫序列{xt}作為輸入變量綠潮覆蓋面積序列{yt}的單輸出傳遞函數(shù)模型：

式中，Ω(B)、E(B)、Θ(B)和Ф(B)為滯后算子B的多項式，其階數(shù)依次為s、r、q和p；b為延遲參數(shù)，即xt的b期滯后值才開始對yt產(chǎn)生影響；εt為隨機干擾性，滿足εt～iidN(0,σ2)，其與xt相互獨立；為傳遞函數(shù)，可表示為V(B)。

3 結(jié)果與討論

3.1 黃海滸苔爆發(fā)期間覆蓋面積和SST變化規(guī)律

2010—2019 年我國黃海綠潮爆發(fā)期內(nèi)滸苔覆蓋面積的時間序列圖呈現(xiàn)出明顯的以年為周期的季節(jié)性單峰波動規(guī)律（見圖1a）。除2010 年和2015年在第8期（7月上旬）達到峰值外，其余年份均在第6 期和第7 期（6 月中下旬）達到峰值，綠潮最大覆蓋面積范圍為330～905 km2，峰值平均為576 km2，峰值中位數(shù)為654 km2。 圖2a 中自相關(guān)函數(shù)（AutoCorrelation Function，ACF）也顯示了綠潮爆發(fā)的循環(huán)特性，綠潮覆蓋面積ACF序列以12期為一個周期，第24 期和第36 期均有明顯的正峰值，第6 期表現(xiàn)出最大負相關(guān)性，說明綠潮覆蓋面積變化總體規(guī)律為1～6期是由初始值逐漸增長到峰值的階段，6～12 期是由峰值逐漸到消亡的階段，這與圖1a 反映出的綠潮覆蓋面積年變化規(guī)律基本一致。

圖1 時間序列圖Fig.1 Time Series Diagram

圖1b顯示，2010—2019年綠潮爆發(fā)區(qū)域海溫變化呈年周期的季節(jié)性增長，SST 范圍為11.2～29.5 ℃，初始溫度區(qū)間為11.2～15.8 ℃，平均溫度為13.8 ℃。2012 年為初始溫度最低值，當年綠潮峰值規(guī)模為10 a 間最小值，2014 年為初始階段溫度最高值，當年綠潮規(guī)模為10 a 間最大值。根據(jù)ACF 統(tǒng)計結(jié)果（見圖2b），將第6 期設定為每年綠潮爆發(fā)的峰值期，峰值階段海溫區(qū)間為19.1～21.5 ℃，平均溫度為20.7 ℃；最終消亡階段（第12期）海溫區(qū)間為25.4～29.5 ℃，平均溫度為27.5 ℃。統(tǒng)計結(jié)果表明，滸苔從生長到消亡的SST 區(qū)間為11.2～29.5 ℃，與吳洪喜等[7]的實驗室海溫梯度（10～30 ℃）對滸苔生長產(chǎn)生影響的分析結(jié)果基本一致。萌發(fā)溫度（第1期）有7 a 低于文獻中滸苔微觀繁殖體的最適宜萌發(fā)溫度15 ℃[5]，也低于陳月紅[11]采用硝酸還原酶活力（Nitric Acid Reductase Activity，NRA）指標測定的滸苔的NRA平均值15 ℃。滸苔快速增長期（4～6 期）的溫度區(qū)間變化范圍為16.4～21.5 ℃，與白雨等[19]的實證研究結(jié)果較為一致，比Taylor 等[8]和Cui 等[21]的實驗室研究成果（15～20 ℃）整體偏高1.5 ℃，在韓紅賓等[5]通過實驗測定的滸苔孢子/配子適宜萌發(fā)溫度15～25 ℃范圍內(nèi)。綠潮出現(xiàn)的最后位置區(qū)域（第12 期）的SST 平均達到27.5 ℃，與韓紅賓等[5]認為的滸苔孢子/配子放散的溫度為20～30 ℃結(jié)論一致，較白雨等[19]實證研究成果偏高1.5℃，與陳月紅[11]梯度水溫在25 ℃組和30 ℃組時滸苔相對生長率下降明顯的結(jié)論相一致。

由圖3可見，綠潮覆蓋面積對不同滯后期（Lag）的SST 序列的互相關(guān)系數(shù)（Cross-Correlation Function，CCF）存在明顯的正弦曲線波動規(guī)律，CCF 超過藍色虛線表明兩序列具有顯著相關(guān)性。Lag=-2 時存在明顯峰值（CCF=0.48），說明SST 序列領先綠潮覆蓋面積序列2 期；Lag=5 時也存在明顯峰值（CCF=-0.46），說明SST 序列還存在滯后綠潮覆蓋面積序列5 期的規(guī)律，CCF 由正轉(zhuǎn)負意味著初期SST 與綠潮覆蓋面積為同方向變化，5 期后反方向變化趨勢最為明顯。

圖3 綠潮覆蓋面積與海溫互相關(guān)系數(shù)圖Fig.3 Correlation coefficient between green tide coverage area and SST

時間序列圖（見圖1）和自相關(guān)函數(shù)圖（見圖2）均表明兩序列具有明顯的季節(jié)周期（s=12），因此分別采用wold 分解法（見圖4）和季節(jié)差分法（見圖5）消除季節(jié)效應。通過wold 分解法提取確定性周期變化規(guī)律，將原序列每期值分別減去每期的平均值來消除序列的季節(jié)性，新生成序列分別記為{area_st}和{temp_st}。季節(jié)差分法是將每期值減去上一年對應期數(shù)據(jù)的差，新生成序列分別記為{area_dt}和{temp_dt}。

圖2 自相關(guān)函數(shù)和偏自相關(guān)函數(shù)圖Fig.2 Autocorrelation function and partial autocorrelation function

圖4 和圖5 的時間序列結(jié)果表明兩種方法均消除了原序列明顯的季節(jié)波動規(guī)律。ACF 結(jié)果表明wold 分解海溫序列呈現(xiàn)出緩慢下降的拖尾趨勢（見圖4a），偏 自 相 關(guān) 函 數(shù)（Partical AutoCorrelation Function，PACF）一階截尾；綠潮覆蓋面積的ACF 和PACF 均呈現(xiàn)出截尾特征（見圖4b）。季節(jié)差分海溫序列也呈現(xiàn)出拖尾特征（見圖5a），并且仍具有正弦波動規(guī)律；綠潮覆蓋面積的ACF 和PACF 均表現(xiàn)出在1 階和12 階顯著（見圖5b），呈一定的季節(jié)波動特征?？傮w來看，wold 分解法去季節(jié)周期效果優(yōu)于差分法，但從ACF 衰減速度來看季節(jié)差分后的序列更為平穩(wěn)。

圖4 wold分解法Fig.4 Wold decomposition method

圖5 季節(jié)差分法Fig.5 Seasonal difference method

3.2 平穩(wěn)性檢驗與最優(yōu)滯后階數(shù)選取

為避免因數(shù)據(jù)出現(xiàn)“偽回歸”而影響估計結(jié)果的有效性，我們首先需要對數(shù)據(jù)進行單位根平穩(wěn)性檢驗和協(xié)整檢驗。綜合采用ADF 檢驗法和PP（Philips & Perron）檢驗法，零假設均為序列為存在單位根的非平穩(wěn)過程，備擇假設為平穩(wěn)序列，在5%的顯著水平下拒絕存在單位根的原假設。結(jié)果表明（見表1），除海溫原始序列{xt}外，其余序列在ADF 檢驗和PP 檢驗下均為無漂移項且無趨勢的平穩(wěn)序列。{xt}序列兩種檢驗結(jié)果均為帶漂移項的平穩(wěn)序列，因此在以{xt}為自變量的回歸方程中引入截距項。由于各序列為I（0）單整序列，原始序列、wold 分解法序列和季節(jié)差分序列之間必然具有長期相關(guān)性的協(xié)整關(guān)系。

表1 單位根檢驗結(jié)果Tab.1 Results of unit root test

自變量和因變量存在協(xié)整關(guān)系，進一步采用Granger 因果檢驗法診斷海溫是否可作為綠潮覆蓋面積模型的自變量，即檢驗海溫是否可作為回歸方程的自變量提高模型預測效果。Granger 檢驗對滯后階數(shù)有較大的敏感性，選取自變量領先因變量0～5 期，原假設為海溫不是綠潮覆蓋面積變化的原因，在5%的顯著水平下拒絕原假設。結(jié)果表明（見表2），原始序列在0～5階均為Granger的統(tǒng)計原因，0階{xt}是{yt}的Granger原因的概率為99.997%，說明海溫當期的變化會引起綠潮覆蓋面積的變化。去掉季節(jié)效應的序列自變量和因變量的Granger 因果檢驗均不能拒絕原假設，說明在模型構(gòu)建時加入海溫數(shù)據(jù)序列不能提高預測精度。

表2 Granger因果關(guān)系檢驗結(jié)果Tab.2 Granger causality test results

3.3 協(xié)整回歸模型的建立和預測

根據(jù)平穩(wěn)性檢驗結(jié)果，{xt}和{yt}均為0 階單整，建立的以{xt}為單變量的回歸模型殘差也具有平穩(wěn)性，符合協(xié)整模型建模要求[20]。采用EG兩步法構(gòu)建協(xié)整模型，通過最小二乘法構(gòu)造回歸模型為：

經(jīng)過ADF 檢驗，式（8）回歸殘差序列{ut}為平穩(wěn)非白噪聲序列，說明兩個變量間存在協(xié)整關(guān)系，但建立協(xié)整回歸模型還需要進一步提取殘差序列中的相關(guān)信息，因此，構(gòu)建殘差序列的ARMA 模型。根據(jù){ut}的ACF 和PACF 的性質(zhì)（見圖6a），回歸殘差序列的ACF 為拖尾特征，PACF 為1 階截尾特征，因此初步設定式（8）中的殘差序列擬合模型為AR（1）。

采用最大似然法對模型參數(shù)進行估計，Box-Ljung 法判斷模型殘差為平穩(wěn)非白噪聲序列。重新進行模型定階數(shù)、參數(shù)估計和殘差檢驗后，構(gòu)建殘差模型階數(shù)為ARMA（2，0，2），各參數(shù)C 均通過顯著性檢驗（見表3），擬合后最終模型殘差為平穩(wěn)白噪聲序列。模型可表示為：

表3 協(xié)整模型殘差參數(shù)顯著性檢驗結(jié)果Tab.3 Significance test results of residual parameters of cointegration model

采用赤池信息準則（AIC）、校正的赤池信息準則（AICc）和貝葉斯準則（BIC）評價模型擬合精度（見表4）。協(xié)整模型擬合精度AIC=13.03，AICc=13.04，BIC=13.22。對2018—2019 年綠潮覆蓋面積進行模型預測（見圖6b），測試集為{yt}序列2018—2019年數(shù)據(jù)。圖6b中黑色線為實際值（2010—2019年），紅色線為擬合值（2010—2017 年）和預測值（2018—2019 年）。采用均方根誤差（Root Mean Square Error，RMSE）和 平 均 絕 對 誤 差（Mean Absolute Error，MAE）評價模型的預測精度（見表4）。協(xié)整模型的預測精度為MAE=122.48，RMSE=171.40。式（9）的計算結(jié)果表明，xt與yt相關(guān)性較強，海溫與綠潮面積變化具有顯著的正相關(guān)性。雖然在6期后綠潮面積有從峰值下降的趨勢，但是6～12期海溫（19.7～29.5 ℃）仍基本處于滸苔最佳生長范圍內(nèi)（10～30 ℃）[7]，僅第12 期平均溫度大于25 ℃，因此式（8）中自變量系數(shù)表現(xiàn)為海溫每升高1 ℃，綠潮覆蓋面積增長15.07 km2。殘差項表明除受海溫影響外，綠潮覆蓋面積還具有顯著的自身變化規(guī)律，滯后因子B 表明綠潮規(guī)模短期變化是ARMA（2，0，2）的自回歸平均移動過程，殘差模型的系數(shù)表明{yt}t-1期對t期的影響大于t-2期的動態(tài)過程，t-1 期的AR 過程和t-2 期的MA 過程參數(shù)為負表明具有負向調(diào)節(jié)趨于平穩(wěn)的作用，6～12 期由于t-1和t-2 期的{yt}較大，負向調(diào)節(jié)功能更為明顯，表現(xiàn)為{yt}的下降趨勢，即綠潮覆蓋面積由峰值衰減的過程。

圖6 協(xié)整模型Fig.6 Cointegration model

表4 協(xié)整模型和傳遞函數(shù)模型精度Tab.4 The precision of co-integration model and transfer function model

3.4 傳遞函數(shù)模型的建立和預測

以{xt}為輸入變量，{yt}為輸出變量構(gòu)建動態(tài)回歸過程的傳遞函數(shù)模型。首先對序列進行預白化處理，基于Box-Jenkins 建模方法對{xt}進行ARMA模型擬合，結(jié)果為帶漂移項的ARMA（2，0，1）模型，傳遞函數(shù)模型參數(shù)通過顯著性檢驗（見表5）和殘差進行白噪聲檢驗，預白化濾波器為：

表5 預白化濾波器參數(shù)顯著性檢驗結(jié)果Tab.5 Significance test results of pre-whitening filter parameters

根據(jù)預白化生成的{xt}殘差序列{αt}和{yt}的派生序列{βt}計算互相關(guān)函數(shù)，排除自變量自相關(guān)干擾后兩者仍具有強相關(guān)性，該結(jié)果與協(xié)整檢驗和因果檢驗結(jié)果保持一致，說明海溫與綠潮覆蓋面積變化的長期趨勢具有強相關(guān)性。由于預白化過程去除了海溫序列的自相關(guān)性，與圖3 的滯后階數(shù)不同，{αt}序列滯后2 期時達到最大值，式（10）中xt的滯后算子為0、B和B2，也表明滯后0～2期海溫對綠潮覆蓋面積有顯著影響。將新生成的{αt}和{βt}進行回歸，模型殘差是平穩(wěn)非白噪聲序列，根據(jù)圖7a將傳遞函數(shù)模型殘差設定為ARMA（0，0，2），整個傳遞函數(shù)模型參數(shù)通過顯著性檢驗（見表6），傳遞函數(shù)模型殘差為平穩(wěn)白噪聲，傳遞函數(shù)模型函數(shù)可表示為：

表6 傳遞函數(shù)模型參數(shù)顯著性檢驗結(jié)果Tab.6 Significance test results of transfer function model parameters

模型擬合精度為AIC=5.40，AICc=5.41，BIC=5.53。對2018—2019 年的綠潮覆蓋面積進行預測（見圖7b），預測精度為RMSE=160.94，MAE=109.70。式（11）中xt的參數(shù)項表明，海溫對綠潮覆蓋面積的影響是顯著的，表現(xiàn)為海溫滯后的0～2期對綠潮覆蓋面積規(guī)模的非線性動態(tài)關(guān)系，其中t 和t-1 期的MA 過程影響最為顯著，t-2 期的AR 過程具有負向調(diào)節(jié)作用。誤差項表明因變量yt除與自變量xt具有相關(guān)性外，還具有綠潮覆蓋面積變化的自身時序關(guān)系，表現(xiàn)為在1～2 期的非線性平均移動回歸過程中，增長階段的yt殘差的負向調(diào)節(jié)機制不明顯，而峰值后t-1 期和t-2 期yt的負向調(diào)節(jié)更為顯著，該結(jié)論與協(xié)整模型中yt的變化規(guī)律一致。

3.5 模型預測精準度評價

將海溫序列進行自回歸預白化處理，提升了模型的擬合精度，即便是對模型復雜性更具有懲罰性的BIC 準則，傳遞函數(shù)模型也優(yōu)于協(xié)整模型。圖6b和圖7b 的擬合結(jié)果表明，從長期趨勢來看，兩個模型均模擬出2013年、2014年和2016年綠潮災害較為嚴重，2012年和2017年綠潮災害較輕。從年內(nèi)模擬效果看，10 a間生消過程估計趨勢與遙感反演結(jié)果趨勢一致，傳遞函數(shù)模型在每年的生消過程中還可以模擬出更多的綠潮覆蓋面積增減的信息。

傳遞函數(shù)模型預測的RMSE 與協(xié)整模型相比提高了6.1%，MAE 提高了10.4%。圖6b 和圖7b 的預測結(jié)果表明，兩模型均預測出了2018年綠潮具有較大規(guī)模，2019 年綠潮規(guī)模處于平均水平，2019 年的傳遞函數(shù)模型預測與真實值非常接近，說明模型具有良好的預測能力。

圖7 傳遞函數(shù)模型Fig.7 Transfer function model

4 結(jié)論和展望

本文基于黃海綠潮遙感覆蓋面積和表層海溫數(shù)據(jù)開展ADF 檢驗和PP 檢驗，結(jié)果顯示兩變量均為0 階單整非白噪聲序列，兩變量間具有長期協(xié)整性。綠潮覆蓋面積自相關(guān)系數(shù)變化規(guī)律表明面積峰值集中在6 月中下旬出現(xiàn)，海溫與綠潮覆蓋的相關(guān)系數(shù)波動規(guī)律表明兩者具有先正相關(guān)后負相關(guān)的季節(jié)特征。分別采用wold 分解法和季節(jié)差分法消除序列中年度波動因素，進一步采用Granger 因果檢驗法發(fā)現(xiàn)，原始海溫序列在領先0～5階均具有顯著互相關(guān)性，說明海溫可作為綠潮覆蓋面積變化模型的自變量因子，消除季節(jié)波動的兩類序列間存在長期協(xié)整性但海溫不能作為綠潮覆蓋面積變化模型的自變量因子。利用協(xié)整模型和傳遞函數(shù)模型，采用海溫因子構(gòu)建靜態(tài)回歸模型和動態(tài)回歸模型，模型參數(shù)的顯著性和殘差檢驗說明兩模型均具有可預測性。綜合采用赤池信息準則和貝葉斯準則評價擬合效果，結(jié)果表明傳遞函數(shù)模型整體優(yōu)于協(xié)整模型。根據(jù)均方根誤差和平均絕對誤差比較模型預測效果，傳遞函數(shù)模型預測精度比協(xié)整模型的RMSE 提升了6.1%，MAE 提升了10.4%。此外，傳遞函數(shù)參數(shù)表明可以利用本周及滯后1～2 周的海溫數(shù)據(jù)預測本周的綠潮覆蓋面積變化，增強了預測的時效性，因此傳遞函數(shù)模型是海溫和綠潮生消變化更適合的預測模型。

綠潮災害涉及多種環(huán)境因子共同作用，其爆發(fā)機制尚未十分清晰，當前從海洋物理、化學和生態(tài)層面預測滸苔綠潮的生態(tài)擴展和消亡過程尚有難度。時間序列方法不要求機制明確，可直接基于歷史數(shù)據(jù)進行綠潮生消模擬和預測，還具有計算量低和業(yè)務性強的特點。本文以海溫因子對滸苔生消變化的影響為研究切入點，得到基于主導因子的綠潮規(guī)模協(xié)整模型和傳遞函數(shù)預測方法，可作為綠潮年度預測的有效途徑。為提高模型預測精準度，下一步還需將光照、營養(yǎng)鹽和降水等其他生態(tài)因子納入模型，并將時間序列模型嵌入漂移擴散預報系統(tǒng)，實現(xiàn)綠潮漂移趨勢和規(guī)模變動預測，以應對防災減災的迫切需要。

致謝：感謝國家衛(wèi)星海洋應用中心為本研究提供了寶貴的基礎觀測資料。