時間序列分析模型在黃海南部小黃魚資源量預測中的應用

宋大德，汪金濤，陳新軍，仲霞銘，熊瑛*，湯建華，吳磊

( 1. 上海海洋大學 海洋科學學院，上海 201306；2. 江蘇省海洋水產研究所，江蘇 南通 226007)

1 引言

小黃魚（Larimichthys polyactis）屬硬骨魚綱（OSTEICHTHYES）、鱸形目（Perciformes）、石首魚科（Sciaenidae）、黃魚屬（Larimichthys），廣泛分布于東海、黃海和渤海近底層水域，也是中、日、韓三國共同利用的主要經濟魚種，在我國近海海洋漁業中具有重要的地位，曾為我國“四大漁業”之一[1]。20 世紀50-80 年代，小黃魚資源經歷了從成熟期到衰退期的演變，主要由于70 年代后捕撈強度增大和環境惡化的影響。20 世紀90 年代后，我國實施了禁漁期、禁漁區和保護區等系列手段以保護小黃魚的產卵場和補充群體，其資源得到了有效保護[2]，產量逐年上升。21 世紀以來，小黃魚產量處于居高不下的態勢，與之相伴隨的是小黃魚低齡化、小型化和性成熟提前等日益凸顯的特征[3–4]，如若不及時控制捕撈強度改變此種狀況，小黃魚漁業很有可能面臨崩潰的風險[5–6]。因此，基于長周期的資源監測調查對小黃魚資源動態進行研究，明確其資源量以有效控制捕撈強度，對當前漁業管理具有重要的參考意義。

迄今為止，在小黃魚資源評估方面缺乏較為系統的定量研究[7–11]，且其中大多是采用基于調查設計的方法，用掃海面積法來評估一個年度春、夏、秋、冬四季資源的時空變化；基于模型的評估方法則較少應用，僅見于剩余產量模型、貝葉斯模型等方法[5,12–13]。而在基于模型的方法中，時間序列分析模型相對于靜態處理統計方法在分析和處理方面具有無可比擬的優勢，隨著數值方法和計算機的發展，已經形成了一套完整的分析和預測體系[14]，在漁情預報漁獲量預測方面已有較好的應用[15–19]，且與人工神經網絡、貝葉斯模型相比，時間序列分析模型表現出更高的精度[20]。

帆式張網（學名：單錨張綱張網）是黃海南部最主要的捕撈形式之一，其作業原理是根據捕撈對象的生活習性和作業區域的水文條件，將網敷設在海中魚、蝦類等的洄游通道上，依靠水流的作用，迫使捕撈對象進入網囊中；其作業位置并非一直固定不變，而是隨著魚群分布變化而轉移作業場所。2008-2015 年《中國漁業統計年鑒》顯示，張網類漁具捕撈量僅次于拖網類、刺網類，位居第三位，占我國近海海洋捕撈量的12%～13%；曾有學者用張網漁獲監測數據估算東海區小黃魚總允許漁獲量[21]和最大持續產量[12]及分析黃海南部小黃魚時空分布特征[2]，以上研究均已為小黃魚漁業資源管理提供了積極的指導作用。本文利用2003-2016 年黃海南部帆式張網小黃魚漁獲監測數據，應用時間序列ARIMA（Autoregressive Integrated Moving Average ，ARIMA）模型，對黃海南部小黃魚資源量進行擬合和預測，旨在為小黃魚的有效管理和可持續利用提供科學依據。

2 材料和方法

2.1 數據來源

數據來源于江蘇省海洋水產研究所2003-2016年1-12 月（伏休期除外）黃海南部帆式張網小黃魚漁獲監測資料。監測船為“蘇啟漁02420”“蘇啟漁03203”“ 蘇 啟 漁02105”“ 蘇 啟 漁02377”“ 蘇 啟漁03650”；所用網具為帆式張網，漁具規格220 m×200 m（網口裝扎周長×縱向網衣拉直長度），最大網目尺寸700 mm（網口），向后沿縱向依次遞減至25 mm（網囊）；2003-2016 年黃海南部帆式張網監測船數、航次數和總投網次數詳見表1。監測范圍涉及黃海南部呂泗漁場、大沙漁場、長江口漁場以及沙外漁場和江外漁場緊鄰125°E 的我國漁船可生產海域（圖1）。

2.2 數據分析

2.2.1 數據處理

單位捕撈努力量漁獲量（Catch Per Unit Effort,CPUE）為每年小黃魚漁獲量除以相應投放網口數，以該CPUE 值作為黃海南部海域小黃魚資源的豐度指標[22]。主要運用SPSS24、Eviews8.0 和Arcmap10.4 軟件進行數據處理。

采用ARIMA 模型對小黃魚年CPUE 進行時間序列擬合及預測，該方法的基本思路是：對于非平穩的時間序列，用若干次階差分使其成為平穩序列，再用ARMA（p，q）對該序列建模。

表1 2003-2016 年黃海南部帆式張網監測船數、航次數和總投網次數Table 1 The number of monitoring vessel by canvas stow net,voyage and hauls from 2003 to 2016 in the southern Yellow Sea

圖1 2003-2016 年黃海南部帆式張網調查區域所涉漁場Fig. 1 The fishing grounds surveyed by canvas stow nets in the southern Yellow Sea from 2003 to 2016

2.2.2 ARIMA（p，d，q）模型

一個ARIMA（p，d，q）過程包括自回歸項（p）、差分階數（d）和移動平均項（q），一般情況下，p、d、q各參數的值不超過2，否則會導致數據失真，用數學公式表示如下：

式中，y(t)、e(t)分別表示原序列和白噪聲序列；B是滯后算子，滿足表達式：Bny(t) =y(t-n)，n= 1，2，···，n；▽d= (1 -B)d是d階差分，d= 1，進行一次差分處理，即令z1(t) = ▽y(t) =y(t) -y(t- 1)；d= 2，進行兩次差分處理，即z2(t) = ▽2y(t) =▽z1(t) =z(t) -z(t-1)，依此類推。

Box-Jenkins 法應用于時間序列建模，預測的過程如下5 個步驟逐級進行：

（1）數據的預處理：包括差分平穩化處理和零均值化處理，這一步驟是使該序列滿足時間序列建模的前提條件。

（2）模型結構的辨識：根據不同模型自相關函數（Autocorrelation Function，ACF）和偏相關函數（Partial Autocorrelation Function，PACF）的性質（表2），確定模型的階數。AR（p）和MA（q）模型的區別在于其所依據的變量不同，AR 模型依據滯后序列進行預測，而MA 模型依據誤差數據序列的未來變化進行預測；多數情況下這兩種模型需要結合使用，就形成了ARMA（p，q）模型。上述3 種模型都是基于平穩的時間序列建立的，當時間序列數據不平穩時，就需要用到ARIMA（p，d，q）模型對數據預處理，如步驟（1）。表2 中，拖尾是指ACF 或PACF 隨著時滯k的增加呈現指數衰減并趨于0，而截尾則是指ACF 或PACF 在某步之后全部為0。

表2 4 類模型的相關函數性質Table 2 The correlation function properties of four types of models

（3）模型參數估計：本文采用最小二乘法估計模型的參數，最小二乘法就是根據一組X和Y的觀測值，選擇Y的逼近函數f（xi；b，c，···）使得偏差的平方和最小，確定逼近函數中的參數b，c，···，以確定Y和X之間的經驗關系。如下式：

式中，b和c1，c2，···，ci代表常數。

（4）模型檢驗：通過對原時間序列與所建模型之間誤差序列是否有隨機性的檢驗來實現；若模型檢驗不能通過，則需要回到步驟（2）重新進行模型結構的鑒別。

（5）模型預測：利用所建立的合適模型導出其預測模型，應用于實際預測。

3 結果

3.1 時間序列的平穩性檢驗

對水平序列、一階差分序列和二階差分序列進行單位根檢驗，結果如表3。

表3 水平序列和差分序列單位根檢驗表Table 3 Unit root test of horizontal sequence and difference sequence

從表3 中可以看出，水平序列的單位根檢驗t值為-2.766 2 明顯大于其1%水平值，水平序列存在單位根，所以水平序列不平穩。經過一階差分后的序列t值略微小于其1%水平值。經過二階差分后的序列t值明顯小于其1%水平值。初步判定模型為ARIMA（p，2，q）模型。

3.2 時間序列的ACF 和PACF 圖

運用SPSS 軟件，繪制出二階差分序列的ACF 和PACF 圖（圖2），由圖2 中ACF 圖可以看出：當延遲期數大于1，其自相關系數位于置信度區間內，因此黃海南部2003-2014 年小黃魚年CPUE 二階差分序列的自相關系數是1 階截尾的。由PACF 圖可以看出其系數值都處在置信度區間內，處于拖尾狀態，因此ARIMA 模型的q值為0。對本研究序列建立ARIMA（1，2，0）和ARIMA（2，2，0）。

3.3 ARIMA 模型擬合及檢驗

利用SPSS 對模型ARIMA（1，2，0）和ARIMA（2，2，0）進行擬合，運行結果顯示：ARIMA（1，2，0）模型貝葉斯信息準則（Bayesian Information Criterion，BIC）值為6.031，ARIMA（2，2，0）模型BIC 值為6.345。依據貝葉斯信息準則，ARIMA（1，2，0）模型擬合的效果較好。利用最小二乘法法對ARIMA（1，2，0）模型進行參數估計（表4）。由表3 得其數學表達式為：yt=-0.714 - 1.006yt-1+et。

通過觀察模型的殘差ACF、PACF 圖（圖3）檢測模型適應性，從圖3 中分析可知，殘差自相關函數和偏相關函數均落在置信區間內，殘差序列為白噪聲過程。因此，用該模型擬合得較好，可將原時間序列中所蘊含的相關性充分提取出來。綜上分析，可確定ARIMA（1，2，0）為最佳模型。

圖2 二階差分序列的自相關函數（ACF）和偏相關函數（PACF）圖Fig. 2 The ACF and PACF diagrams of second order difference sequence

表4 ARIMA（1，2，0）模型參數估計表Table 4 The parameter estimation of the ARIMA (1, 2, 0) model

圖3 ARIMA（1，2，0）模型殘差序列自相關函數（ACF）和偏相關函數（PACF）圖Fig. 3 The ACF and PACF diagrams of the ARIMA (1, 2, 0)model residual sequence

3.4 ARIMA 模型預測

運用ARIMA（1，2，0）模型，對2003-2014 年黃海南部小黃魚年CPUE 時間序列進行擬合，結果如圖4，擬合值與實際值相關系數為0.881，且相關性顯著（p＜0.05）。對2015-2016 年黃海南部小黃魚年CPUE值進行預測，模型預測值分別為47.66 kg/網和49.16 kg/網，與實際值51.10 kg/網和40.05 kg/網的相對誤差分別為6.73% 和22.75%，總體相對誤差為14.74%。因此，可說明運用ARIMA（1，2，0）模型能有效地對黃海南部小黃魚年CPUE 擬合和預測。

圖4 2003-2014 年黃海南部小黃魚年CPUE 模型擬合值與實際值關系Fig. 4 The relationship between simulated and actual CPUE of small yellow croaker in the southern Yellow Sea from 2003 to 2014

4 分析與討論

2003-2016 年黃海南部帆式張網的監測數據顯示，小黃魚年CPUE 值的變化趨于平穩，表明近14 年小黃魚漁獲產量并沒有出現下降的態勢。這與《中國漁業統計年鑒》中江蘇海域小黃魚產量一直居高不下的趨勢一致，亦與嚴利平等[23]在小黃魚資源研究中認為“小黃魚資源自20 世紀90 年代后半期以來處于穩定豐厚期”的結論一致，究其原因：其一，高強度捕撈是造成黃海小黃魚漁獲量較高的主要因素[24]；其二，我國政府自1995 年以來實施伏季休漁制度并不斷調整優化，并于2010 年設立“呂泗漁場小黃魚銀鯧國家級水產資源保護區”等舉措對小黃魚資源起到增殖效果[2,23]；其三，張吉喆[25]認為中韓兩國2001 年簽署的《中華人民共和國政府和大韓民國政府漁業協定》，對小黃魚資源也起到了很好的養護作用；其四，小黃魚1 齡魚性成熟比例和個體絕對繁殖力顯著提高，這對小黃魚的資源補充能力提升起到了關鍵作用[6]。過度捕撈是誘導小黃魚生命史特征演變的關鍵因素，許多開發性魚類資源因過度捕撈而衰退甚至枯竭[26]。鑒于此，應加強小黃魚資源評估，選取最優模型和參數配置對管理措施下的現行資源狀況進行正確評估和精準預測，以確保真正實現漁業資源種群結構的好轉和合理化。

通過比較掃海面積法、剩余產量模型、貝葉斯模型和時間序列模型在小黃魚資源評估上擬合和預測效果，逐一分析了各類評估方法和模型使用時的缺陷，以及影響評估效果的因素。采用掃海面積法評估小黃魚資源狀態的文章較少，林龍山[8]曾用掃海面積法對東海區小黃魚資源量進行估算，但實際產量與估算結果相差較大；關麗莎等[27]曾用掃海面積法和地階廣義線性混合模型分別對黃海南部小黃魚資源量進行評估，結果表明模型估算精度比掃海面積法估算的要高。分析其原因認為，掃海面積法中所使用小黃魚的可捕系數不僅受其晝夜垂直移動特征影響[27]，還與棲息水域、首次性成熟年齡[28]以及個體大小[29]等有關。

多數研究者采用剩余產量模型對小黃魚資源狀況進行評估，其中參數估計是影響模型預測精度的主要因素，如自然死亡系數M 對生長方程中的參數、生物參考點估計影響較大，進而造成模型預測結果的偏差[30]。林龍山[8]研究2002 年東海區域小黃魚資源，估算M 值為0.58；周永東等[12]研究1996-2009 年東海區域小黃魚資源，估算M 值約為0.12；高春霞等[31]研究2016 年浙江南部近海小黃魚資源時估算M 值為1.34。上述研究中小黃魚自然死亡系數M 值偏離較大，筆者認為，除調查水域小黃魚棲息水溫因素影響外，更有可能是估算M 值的方程由于調查區域不同、年代變遷等因素出現不適用性。例如，詹秉義等[32]曾在1984 年評估綠鰭馬面鲀資源時推導出M 值的線性回歸方程，該方程在多種漁業資源群體中一直沿用至今，其中包括小黃魚群體?；谏鲜鯩 值變化的分析，筆者認為該線性方程可能已不適用于當前我國小黃魚資源群體評估。日本學者田中昌一提出采用年齡數據估算魚類自然死亡系數 M 值的線性回歸方程需要更多的、更新的數據進行及時修正，否則參數設置所得出的評估結果作為漁業管理參考依據將會影響漁業資源可持續利用管理決策的正確性[33]。官文江和吳佳文[34]采用剩余產量模型對印度洋黃鰭金槍魚資源評估時發現，形狀參數很難被估計，而且隨著數據時間段變化，形狀參數的最佳取值差別很大。使用模型過程中參數估算的不確定性是無法避免的，需采用更為科學準確的方法降低不確定性，此后有學者開始采用貝葉斯方法試圖解決這一問題[35–37]，貝葉斯方法充分考慮了模型和參數值存在的不確定性，其參數后驗分布集合了其他參數所有可能值，因此可以提高模型評估的可信度。

時間序列是一組依賴于時間t 的隨機變量，這組隨機變量所具有的自相關性和依存關系表征了預測對象發展的延續性，用相應的數學模型近似描述這種相關性，就可以從時間序列的過去值，預測其將來的值[38]。時間序列分析是一種數據驅動的方法，Bakker 等[39]采用多種模型研究地下水流問題時，發現時間序列分析模型較常規評估模型不僅簡單，而且擬合效果更好。董江水和諸英富[16]曾基于統計數據用時間序列分析模型對江蘇省河蟹總產量進行預測，模型擬合值與實際值的相關系數為0.995，產量預測相對誤差為3.8%，相較本文模型擬合與預測結果要高，其原因可能是本文的原始時間序列數據取自固定監測樣本船，相較統計數據波動性較大，因此模型的擬合值與實際值間的相關性就會被削弱，致使預測精度下降。但本模型預測精度近90%，進一步表明帆式張網監測船漁獲數據適用于時間序列分析模型的預測。

需要注意的是，即使是同一海域、同一物種，同一時間序列數據也不能簡單套用時間序列分析模型，不同時間序列數據ARIMA 模型的p、d、q 值不盡一致，應根據相關理論指導和分析，對具體的資料進行詳盡分析以確定適宜的p、d、q 值[40]。本文構建的ARIMA（1，2，0）模型通過了適應性檢驗，能較好地擬合黃海南部小黃魚資源量的變動趨勢，說明時間序列模型用于漁業資源量預測是可行的，但由于ARIMA模型的局限性，在短期預測方面更具優勢[38]。