西北太平洋柔魚豐度的灰色災變預測

解明陽，陳新軍,2,3,4,5*

(1. 上海海洋大學 海洋科學學院，上海 201306；2. 農業(yè)農村部大洋漁業(yè)開發(fā)重點實驗室，上海 201306；3. 國家遠洋漁業(yè)工程技術研究中心，上海 201306；4. 大洋漁業(yè)資源可持續(xù)開發(fā)教育部重點實驗室，上海 201306；5. 農業(yè)農村部大洋漁業(yè)資源環(huán)境科學觀測實驗站，上海 201306)

1 引言

柔魚（Ommastrephes bartramii）是大洋經濟頭足類，廣泛分布于西北太平洋，是我國魷釣漁船的重要捕撈對象[1]。自1993年我國開始開發(fā)柔魚資源以來，捕撈規(guī)模與作業(yè)海域不斷擴大，年產量基本穩(wěn)定在6萬～10萬t[2]。由于柔魚是一種“生態(tài)機會主義”物種，其中心漁場與資源豐度極易受海洋環(huán)境與氣候變化的影響[3]，因此柔魚資源量會出現劇烈的年間變化。已有研究表明，柔魚的資源量變化與環(huán)境因子變化，例如海表溫度（Sea Surface Temperature, SST）[4]、黑潮[5]有密切的關系；資源量的驟增驟減現象也受氣候變化的影響，例如 El Ni?o-La Ni?a 現象[6]、太平洋年代際濤動（Pacific Decadal Oscillation，PDO）現象[7]等影響。目前國內外學者已經對北太平洋柔魚漁場時空分布[8]、資源評估與管理[9]、資源量預報模型[10]等做了較為全面的研究。對柔魚資源量進行年度的災變預測，將有利于掌握西北太平洋柔魚的資源狀況，為資源合理利用提供依據。

灰色系統(tǒng)理論是一門有關不確定系統(tǒng)的理論，其優(yōu)點在于允許樣本數量較少且服從任意分布，這對于缺乏數據的漁業(yè)科學領域來說有很好的應用前景[11]。目前，灰色系統(tǒng)理論在資源量預測[12]、漁場空間分布[13]、漁業(yè)水質評價[14]等方面已取得了較好的成果。為此，本研究采用灰色系統(tǒng)研究方法中的灰色災變預測模型，對西北太平洋柔魚資源豐度的豐歉年進行預測，以達到預警目的，為北太平洋漁業(yè)生產企業(yè)和管理部門及時進行有效管理提供技術支撐。

2 材料與方法

2.1 漁業(yè)生產數據

漁業(yè)生產統(tǒng)計數據來自上海海洋大學中國魷釣技術組，時間為1995-2017年，空間范圍為35°～45°N、140°～179°E，統(tǒng)計內容包括日期、經度、緯度、日產量和作業(yè)漁船，空間分辨率為1°×1°。以單位捕撈努力量漁獲量（Catch Per Unit Effort，CPUE）表征柔魚資源豐度，以1°×1°為一個漁區(qū)，計算各年每個漁區(qū)的CPUE值。

SST數據來源于美國國家海洋和大氣管理局（National Oceanic and Atmospheric Administration,NOAA）的網站（http://oceanwatch.pifsc.noaa.gov/erddap/index.html），時間范圍為 1995-2017 年的1-12月，時間分辨率為月，利用平均法計算各年作業(yè)漁場平均SST。

2.2 CPUE標準化

統(tǒng)計每年所得到的總漁獲量（Catch），以每年捕撈的船數作為捕撈努力量(Effort)，計算名義CPUE，其公式為

式中，CPUE為每年的單位捕撈努力量，單位為t/（船·a）；Catch的單位為 t/a。

利用廣義線性回歸模型（Generalized Linear Model，GLM）對名義CPUE進行標準化。由于CPUE受環(huán)境因子（SST）、漁場空間位置（漁場重心的經度、緯度）等因素影響，所以將其作為解釋變量，CPUE作為響應變量，去除空間和環(huán)境因子對CPUE的影響，以名義CPUE和標準化的CPUE作為柔魚資源豐度的指標。

漁場重心計算公式為

式中，Lon、Lat分別為漁場重心位置的經度和緯度；Ci為i漁區(qū)的產量；Loni、Lati分別為漁區(qū)i的中心經緯度位置；K為漁區(qū)的總個數。

GLM的表達形式為：glm[CPUE]～SST +Lon+Lat+ε。其中，SST為年平均海表面溫度，ε為誤差項。

2.3 模型的構建

利用GM（1, 1）模型對北太平洋柔魚資源豐度進行災變預測[15]。具體方法如下：

設原始序列X=（x(1)，x(2)，···，x(n)），取序列的第3分位數Q3（大于75%的值）為上限災變值 ξ1，取序列第1分位數Q1（小于25%的值）為下限災變值 ξ2，將序列X中大于或小于 ξ的點認為是異常值的點[16]，將其組成上限或下限災變序列Xξ=（x[q(1)]，x[q(2)]，···，x[q(m)]），設災變日期序列Q(0)=（q(1)，q(2)，···，q(m)），對其建立GM（1, 1）模型進行災變預測。GM（1, 1）模型計算方法如下[17]：

令Q(1)=（x(1)(1)，x(1)(2)，···，x(1)(n)）為災變日期序列Q(0)的一次累加生成序列。其中令Z(1)=（Z(1)(2)，Z(1)(3)，···，Z(1)(n)）為Q(1)的緊鄰均值生成數列，其中，k=2,3,···,n，得到GM（1, 1）模型的基本形式為：Q(0)(k)+aZ(1)(k)=b，a為發(fā)展系數。

2.4 模型的檢驗

灰色災變GM（1, 1）模型的驗證主要包括殘差檢驗與后驗差檢驗[18]。殘差為原始值與預測值之間的絕對誤差與相對誤差；后驗差檢驗是按照殘差的概率分布進行檢驗，記原始序列標準差絕對誤差序列標準差方差比小誤差概率

殘差越小，模型精度越好。方差比C越小，說明S2小，S1大，即殘差方差小，原始數據方差大，殘差比較集中，擺動幅度小，原始數據比較分散，擺動幅度大，所以模型模擬效果更好。模型精度檢驗等級見表1。

發(fā)展系數a反映主因素與各子因素之間的協調程度，當-a≤0.3時，模型可用于中長期預測；當0.3＜-a≤0.5時可用于短期預測[19]。

表1 灰色系統(tǒng)理論模型精度檢驗[18]Table 1 Precision test accuracy of grey model[18]

3 研究結果

3.1 CPUE情況與災變點

1995-2017年 各年的名義平均CPUE與標準化CPUE見圖1。經過GLM模型標準化的CPUE與名義CPUE的變化趨勢基本一致但波動較小。2007年名義 CPUE 值最大，為 4.197 t/（船·a），2009 年名義CPUE 最小，為 1.319 t/（船·a）；2017 年標準化的CPUE 值最大，為 2.928 t/（船·a），2002 年標準化的CPUE 最小，為 1.446 t/（船·a）。

根據災變值的劃分（下限災變值為2.13 t/（船·a），上限災變值為 2.39 t/（船·a）），以 GLM 模型標準化CPUE建立的災變模型的上限災變點為2005年、2006年、2007年、2008年、2010年、2017年，災變序列為 11、12、13、14、16、23；下限災變點為 1996年、1999年、2001年、2002年、2015年、2016年，災變序列為 2、5、7、8、21、22。

3.2 模型的建立與檢驗

從模型精度來看，以GLM模型標準化CPUE建立的下限災變預測模型的平均相對誤差為15.32%，上限災變預測模型的平均相對誤差為8.19%（表2）。從發(fā)展系數a來看，下限災變預測模型（a=-0.064）和上限災變預測模型（a=0.117）均可用于中長期預測（表3）。從后驗差檢驗結果（表3）來看，下限災變預測模型（C=-0.342，P=1）和上限災變預測模型（C=-0.341，P=1）的方差比C小于 0.35，小誤差概率P大于0.95，模型精度等級為Ⅰ級。

根據表3的響應函數來看，下一次的下限災變和上限災變發(fā)生序號分別為34.86和25.20，在2016年以后的11.14 a（由34.86 a減去23.72 a得出），即2027年會發(fā)生資源豐度歉年；在2017年以后的4.0 a（由25.20 a減去21.20 a得出），即2021年會發(fā)生資源豐度豐年（表 3）。

4 討論與分析

4.1 氣候和環(huán)境變化對西北太平洋柔魚豐度的影響

圖1 1995-2017年西北太平洋柔魚資源豐度變化Fig. 1 CPUE of Ommastrephes bartramii in the Northwest Pacific from 1995 to 2017

表2 灰色災變GM（1, 1）模型預測值與真實值之間的相對誤差Table 2 Relative error between the predicted value and the true value of the GM（1, 1）model

表3 灰色災變GM（1, 1）模型的相關參數Table 3 Related parameters of the GM（1, 1）model

柔魚作為短生命周期的種類，一生產卵一次，因此其漁獲量的多少取決于當年的柔魚補充量[20]。造成柔魚資源豐度下降的原因是多方面的，如沈建華等[21]認為2001年低產的原因主要為：黑潮大蛇行、黑潮勢力偏弱、北部水溫過低、漁汛期流隔不明顯和餌料不充足。陳峰等[22]認為2009年低產的原因主要為產卵場發(fā)生黑潮大彎曲，傳統(tǒng)作業(yè)漁場有一個明顯的冷水南下。唐峰華等[23]認為SST偏低導致的補充量減少，葉綠素a濃度分布與海流波動較大導致的餌料分布不穩(wěn)定等現象也是2001年低產的原因之一。從大尺度氣候變化來看，PDO是一種以10 a周期變化的太平洋氣候變化，PDO指數定義為20°N以北北太平洋SST異常第一主分量。余為等[7]對不同氣候模態(tài)下西北太平洋柔魚漁場環(huán)境做出分析，認為PDO暖期會為仔幼魚提供較好的攝食環(huán)境，易形成適宜的育肥場，導致CPUE較高，PDO冷期則相反。根據美國華盛頓大氣與海洋聯合研究所的數據（http://jisao.washington.edu/pdo/PDO.latest）結合災變年份的分布來看（圖2），基本上災變年份符合PDO冷暖期的分布規(guī)律，即上限災變年份位于PDO暖期，下限災變年份位于PDO冷期。而且在冷暖期轉換的年份（2001-2002年，2007-2008年，2014-2015年）災變年份類型仍維持在轉變期前的狀態(tài)，說明PDO對柔魚資源的影響具有一定的滯后性。此外，El Ni?o-La Ni?a事件是太平洋海域氣候變率的最強信號，對柔魚資源豐度也有較大的影響。如Chen等[6]認為La Ni?a事件發(fā)生時，柔魚產卵場溫度升高，不利于仔幼魚的繁殖與孵化，導致資源補充量下降，反之El Ni?o事件時資源補充量上升。Yu等[24]通過構建棲息地模型發(fā)現La Ni?a事件下產生的暖水團，使適宜的葉綠素濃度和海表面高度距平范圍擴大，導致棲息地面積擴張，而El Ni?o事件則相反。根據NOAA對El Ni?o和 La Ni?a事件的定義（http://ggweather.com/enso/oni.htm），即若Ni?o3.4區(qū)SST異常連續(xù)5個月滑動平均值超過+0.5℃認為發(fā)生一次El Ni?o事件；若連續(xù)5個月低于-0.5℃則認為發(fā)生一次La Ni?a事件。結合災變年份分布來看（圖3），基本上上限災變年份發(fā)生在El Ni?o事件期間，下限災變年份發(fā)生在La Ni?a事件期間。相比于 PDO 而言，El Ni?o 和 La Ni?a事件期間的災變年份分布更為明顯，如下限災變年1996年對應La Ni?a事件，上限災變年2003-2005年對應El Ni?o年，并且其滯后性比PDO更小。

4.2 災變預測模型分析

災變預測年份顯示（表3），在4 a或者更短的時間內，北太平洋柔魚的資源豐度依舊很高，這也符合今后持續(xù)的PDO暖期與El Ni?o事件的氣候變化背景，但在 PDO 冷暖期與 El Ni?o-La Ni?a事件轉化的年份，易發(fā)生資源豐度歉年，并在11 a左右的下一個周期開始仍處于低資源豐度的狀態(tài)。所以建議北太平洋漁業(yè)生產企業(yè)和管理部門在PDO暖期與El Ni?o事件期間加大生產與捕撈努力量的投入，爭取經濟效益的最大化；相反，在PDO冷期與La Ni?a事件期間及時做好低產的準備，減少船只、人力、管理等方面由于多余投入而帶來的經濟損失。

圖2 1995-2017年太平洋年代際濤動（PDO）指數時間序列變化及災變年份分布Fig. 2 Time series of PDO and catastrophe year distribution from 1995 to 2017

圖3 1995-2017年El Ni?o和La Ni?a事件變化及災變年份分布Fig. 3 El Ni?o and La Ni?a event changes and catastrophic years distribution from 1995 to 2017

灰色災變預測實質上異常值預測，而異常值的設置往往是根據人們憑經驗確定的，具有很強的主觀性，因此災變閾值的確定十分重要。本研究采用統(tǒng)計學上的四分位數法確定災變閾值的上下限，從模型的精度來看，災變GM（1, 1）模型可以有效地預測災變發(fā)生時間。從預測模型的結果看（表2），對于短期的預測值精度都較高，但隨著預測時序的增長，會出現誤差較大的點，證明灰色災變預測模型僅適合短期的災變點預測。此外，災變預測不能預測連續(xù)發(fā)生的災變點，如本研究的連續(xù)災變年份包括上限災變點2005-2008年，下限災變點2000-2001年，2015-2016年。由于GM（1, 1）模型過于依賴樣本序列，所以樣本長度不同，序列起始劃分不同，預測值會具有一定的不確定性，如樣本太少會導致結果失真，樣本太多難以達到所要求序列滿足的光滑度。建議以后的研究采用多樣本序列組合形式建立多種模型，考慮漁獲量或捕撈努力量等指標，篩選出最優(yōu)模型用于災變長期預報。