基于環(huán)境因子的東南太平洋智利竹筴魚剩余產(chǎn)量模型建立

白思琦， 鄒曉榮，2*， 丁鵬， 林銘

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

智利竹筴魚（Trachurus murphyi）是高度洄游的中上層海洋經(jīng)濟魚種[1-2]，具有較高商業(yè)價值，主要分布在智利、厄瓜多爾、秘魯?shù)妊睾蛞约爸抢麑俳?jīng)濟區(qū)向西延伸的外海海區(qū)[2-3]。自20世紀70年代以來，其產(chǎn)量多年位居世界單一魚種前列，在一些遠洋漁業(yè)國家和相關(guān)沿海國家的捕撈量中占有較大比重[4-7]。1995年產(chǎn)量達到最高峰，為4.96×106t[8]。進入21世紀后，智利竹筴魚資源量受到氣候的長期波動和人類活動的影響，出現(xiàn)較大震蕩，漁場出現(xiàn)明顯變動[7]。

已有學(xué)者針對海洋環(huán)境變化對智利竹筴魚資源的影響進行了研究。楊香帥等[8]采用灰色關(guān)聯(lián)度分析方法、Li 等[9]通過不同權(quán)重的棲息地適應(yīng)性指數(shù)模型、Laura 等[10]通過神經(jīng)網(wǎng)絡(luò)模型、汪金濤等[11]和陳春光[12]通過主成分和BP 神經(jīng)網(wǎng)絡(luò)方法、晉偉紅[13]通過偏最小二乘法計算模型等以及楊嘉樑等[14]根據(jù)海面溫度、葉綠素a 濃度、海面高度對智利竹筴魚海洋環(huán)境因子與中心漁場變動及資源的關(guān)系進行了研究。但根據(jù)多環(huán)境因子預(yù)測智利竹筴魚產(chǎn)量的研究鮮少報道。李綱等[15]、鄭曉瓊等[16]與徐香香等[17]通過建立基于環(huán)境因子的剩余產(chǎn)量模型（surplus production model）預(yù)測魚種資源量。剩余產(chǎn)量模型對數(shù)據(jù)要求不高且簡單實用，在漁業(yè)中應(yīng)用廣泛，但傳統(tǒng)的剩余產(chǎn)量模型無法將環(huán)境因子應(yīng)用于漁業(yè)資源的預(yù)測。本研究根據(jù)漁撈日志數(shù)據(jù)，基于海洋環(huán)境因子的剩余產(chǎn)量（environmentally dependent surplus production，EDSP）模型，將海面溫度（sea surface temperature，SST）、葉綠素 a 濃度（chlorophyll a concentration，Chl-a）、海面溫度異常值（sea surface temperature anomaly，SSTA）和海面高度異常值（sea surface height anomaly，SSHA）4 種環(huán)境變量引入模型中，通過環(huán)境因子不同的組合方式對智利竹筴魚漁業(yè)資源豐度進行動態(tài)研究和預(yù)測，探索環(huán)境變化對智利竹筴魚資源豐度產(chǎn)生的影響，為智利竹筴魚的合理捕撈提供參考和幫助。

1 材料與方法

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

智利竹筴魚捕撈數(shù)據(jù)來源于中國遠洋漁業(yè)協(xié)會大型加工拖網(wǎng)項目技術(shù)組提供的漁撈日志，捕撈數(shù)據(jù)包括作業(yè)時間、起放網(wǎng)時間、作業(yè)位置以及漁獲量等。作業(yè)區(qū)域為19°—46°S、73°—92°W，時間范圍選擇為2012—2018年4—8月，空間分辨率為0.5°×0.5°。

環(huán)境數(shù)據(jù)SST、Chl-a、SSTA、SSHA 均來源于美 國 國 家 大 氣 局（National Oceanic and Atmospheric Administration，NOAA）Oceanwatch 網(wǎng)站（www.oceanwatch.pifsc.noaa.gov），時間范圍為2012—2018年 4—8月，空間分辨率為 0.125°×0.125°，時間分辨率為周。

1.2 數(shù)據(jù)處理及分析

1.2.1 單位捕撈努力量漁獲量計算 單位捕撈努力量漁獲量（catch per unit effort，CPUE）是指在一定時間內(nèi)，某漁場平均每個單位捕撈努力量所能獲得的漁獲量，是用于表示智利竹筴魚資源豐度的主要指標，其計算公式[7]如下。

式中，Ci表示從第i天的累計捕撈量（t）；Ti表示從第i天的捕撈時間（h）；YCPUE為單位捕撈努力量漁獲量（t·h-1）。

1.2.2 CPUE 標準化 CPUE 標準化是漁業(yè)資源評估必需的指標[18]，目前我國漁業(yè)CPUE 標準化應(yīng)用最為廣泛的模型為廣義加性模型（generalized additive model，GAM）[19-21]，本文使用該模型對智利竹筴魚CPUE 進行標準化，以去除相關(guān)因子對智利竹筴魚CPUE 的影響。將SST、Chl-a、SSTA 和SSHA 4個環(huán)境因子和經(jīng)緯度作為解釋變量加入到GAM 模型中，并對其進行分析，GAM模型表達式[22]如下。

式中，ε為誤差項。在用R 軟件處理時，為避免CPUE 出現(xiàn)零值報錯，故將所有CPUE 結(jié)果加1后再處理。

1.2.3 關(guān)鍵因子選取 利用SPSS 23軟件對2012—2018年智利竹筴魚漁場范圍內(nèi)分辨率0.5°×0.5°的SST、Chl-a、SSTA 和 SSHA 環(huán)境數(shù)據(jù)和標準化后的CPUE 序列值進行皮爾遜相關(guān)性分析，得出與CPUE 相關(guān)性最高的環(huán)境因子，將CPUE 殘差與該環(huán)境因子數(shù)據(jù)進行線性回歸分析，選出關(guān)聯(lián)度最高的海域，并將該海域環(huán)境數(shù)據(jù)作為表征影響智利竹筴魚資源豐度的指標，建立EDSP模型。

1.2.4 基于環(huán)境因子的剩余產(chǎn)量模型 傳統(tǒng)的剩余產(chǎn)量模型一般用于資源評估[23]，將環(huán)境因子引入剩余產(chǎn)量模型中建立EDSP 模型可以用于研究環(huán)境變化對漁業(yè)資源的影響，可以有效地促進漁業(yè)資源科學(xué)研究及合理利用。

本文建立5種EDSP模型，環(huán)境因子分配如下：①只考慮SST對智利竹筴魚資源量的影響；②同時考慮SST 和Chl-a 對智利竹筴魚資源量的影響；③同時考慮SST、Chl-a和SSHA對智利竹筴魚資源量的影響；④同時考慮SST、Chl-a 和SSTA 對智利竹筴魚資源量的影響；⑤同時考慮SST、Chl-a、SSTA和SSHA對智利竹筴魚資源量的影響。一般剩余產(chǎn)量模型計算公式[15]如下。

式中，Cn為第n月的漁獲量；En為第n月的標準捕撈努力量；Un為第n月的標準 CPUE；a、b 為參數(shù)。

假設(shè)智利竹筴魚漁獲量的實際測量值C和預(yù)測估計值c的差為標準捕撈努力量的比值，即CPUE 的月間變化ΔUn是由環(huán)境變化造成的，可建立如下方程。

式中，G(Fn)為第n月由關(guān)鍵環(huán)境因子導(dǎo)致的CPUE波動；φ為常數(shù)。

根據(jù)公式（3）～（5）得到環(huán)境因子與智利竹筴魚漁獲量的關(guān)系模型如下。

式中，Δun為ΔUn的預(yù)測值。

1.2.5 模型擬合度評價 對EDSP 模型誤差項進行顯著性檢驗，常用的方法為赤池信息準則(akaike information criterion，AIC)檢驗[24]、蒙特卡羅檢驗[25]和F檢驗[26-27]等。本文選取 AIC 檢驗，該方法多用于比較不同模型的顯著性，當(dāng)2個模型AIC值相差超過3 時，較低AIC 值的模型擬合效果更好，表達式如下。

式中，n為樣本數(shù)；RSS 為殘差平方和；P為模型中參數(shù)個數(shù)。

本文使用平均絕對誤差（mean absolute error，MAE）和平均相對誤差（mean relative error，MRE）對建立的模型預(yù)測結(jié)果精度進行評價[17]，誤差越小則模型預(yù)測精度越高。

式中，n為樣本數(shù)；ca為實際漁獲量，ce為預(yù)測漁獲量。

1.3 數(shù)據(jù)統(tǒng)計分析

本文采用SPSS 23 軟件對數(shù)據(jù)進行統(tǒng)計分析，使用Excel軟件進行圖表的繪制。

2 結(jié)果與分析

2.1 相關(guān)性分析及關(guān)鍵因子的選取

對2012—2018年智利竹筴魚漁場范圍內(nèi)分辨率0.5°×0.5°的SST、Chl-a、SSTA 和 SSHA 環(huán)境數(shù)據(jù)和CPUE 序列值進行相關(guān)性分析，SST、Chl-a、SSTA 與 CPUE 的相關(guān)系數(shù)分別為 0.035、0.049、0.041，相關(guān)性顯著（P＜0.05），環(huán)境因子SSHA 與CPUE 的相關(guān)性（0.055）較弱，其中相關(guān)性最高的環(huán)境因子為SST，將該環(huán)境數(shù)據(jù)與CPUE 殘差做回歸分析，得到關(guān)聯(lián)度最高的海域為40.5°S、79.0°W（表1），因此本研究將40.5°—41.0°S、79.0—79.5°W 海域作為表征影響智利竹筴魚資源豐度的指標。

表1 CPUE殘差回歸分析Table 1 Regression of the residual of CPUE

2.2 智利竹筴魚一般剩余產(chǎn)量模型

對智利竹筴魚2012—2018年漁獲量以及標準化CPUE 數(shù)據(jù)進行線性回歸分析（圖1），得到剩余產(chǎn)量模型參數(shù)a 和b，數(shù)值分別為1.3083 39 和0.000 13，根據(jù)式（3）（4），可得智利竹筴魚一般剩余產(chǎn)量模型的表達式如下。

圖1 智利竹筴魚CPUE與漁獲量的關(guān)系Fig.1 Relationship between the catch per unit effort of Chilean jack mackerel and catch

2.3 基于環(huán)境因子的智利竹筴魚剩余產(chǎn)量模型分析

2.3.1 基于SST的智利竹筴魚剩余產(chǎn)量模型 SST對于CPUE的影響最為顯著，選取SST為關(guān)鍵環(huán)境因子建立剩余產(chǎn)量模型EDSP1 模型。標準CPUE殘差與各月海域40.5—41.0 °S、79.0—79.5 °W位置對應(yīng)的SST 的擬合曲線如圖2 所示?；貧w分析表明，智利竹筴魚標準CPUE殘差與各月海域位置對應(yīng)的SST呈線性關(guān)系，根據(jù)式（5）可得式（11）。

圖2 智利竹筴魚CPUE殘差與40.5°S、79.0°W海域SST的擬合曲線Fig.2 Fitted quadratic curve of residual of CPUE of Chilean jack mackerel against monthly SST at 40.5°S，79.0°W

當(dāng) ΔUn=0、Tn=14.35 ℃，由圖2 可知，當(dāng)對應(yīng)月份海域40.5—41.0°S、79.0—79.5°W 位置的 SST低于14.35 ℃，智利竹筴魚標準CPUE 殘差小于0；當(dāng)對應(yīng)月份SST 高于14.35 ℃，智利竹筴魚標準CPUE 殘差大于0，最適SST 為14～16 ℃。根據(jù)式（6）可得基于SST 的剩余產(chǎn)量模型EDSP1 的最終表達式。

2.3.2 基于SST、Chl-a 的智利竹筴魚剩余產(chǎn)量模型 顯著性分析結(jié)果表明，環(huán)境因子Chl-a對于智利竹筴魚CPUE 產(chǎn)生一定影響，故關(guān)鍵環(huán)境因子加入Chl-a，建立SST 和Chl-a 的剩余產(chǎn)量模型EDSP2。智利竹筴魚標準CPUE 殘差與各月海域40.5°S、79.0°W 位置對應(yīng)的 SST 和 Chl-a 進行回歸分析。

根據(jù)式（6）可得基于SST 和Chl-a 的剩余產(chǎn)量模型EDSP2的最終表達式。

2.3.3 基于SST、Chl-a 和SSHA 的智利竹筴魚剩余產(chǎn)量模型 SSHA 對于智利竹筴魚CPUE 產(chǎn)生的影響小于Chl-a，故在第3個剩余產(chǎn)量模型中加入關(guān)鍵環(huán)境因子SSHA，建立SST、Chl-a 和SSHA的剩余產(chǎn)量模型EDSP3。智利竹筴魚標準CPUE殘差與各月海域40.5°S、79.0°W 位置對應(yīng)的SST、Chl-a和SSHA進行回歸分析。

根據(jù)式（6）可得基于SST、Chl-a 和SSHA 的剩余產(chǎn)量模型EDSP3的最終表達式。

2.3.4 基于SST、Chl-a和SSTA的智利竹筴魚剩余產(chǎn)量模型 顯著性分析結(jié)果表明環(huán)境因子SSTA對于智利竹筴魚CPUE 產(chǎn)生的影響小于Chl-a，故在第4個剩余產(chǎn)量模型中加入關(guān)鍵環(huán)境因子SSTA，建立 SST、Chl-a 和 SSTA 的剩余產(chǎn)量模型EDSP4。智利竹筴魚標準CPUE 殘差與各月海域40.5—41.0°S、79.0—79.5°W 位置對應(yīng)的SST、Chla和SSTA進行回歸分析.

根據(jù)式（6）可得基于SST、Chl-a 和SSTA 的剩余產(chǎn)量模型EDSP4的最終表達式。

2.3.5 基于環(huán)境因子 SST、Chl-a、SSTA 和 SSHA 的智利竹筴魚剩余產(chǎn)量模型 將4個環(huán)境因子均加入剩余產(chǎn)量模型中，建立SST、Chl-a、SSTA和SSHA的剩余產(chǎn)量模型EDSP5。智利竹筴魚標準CPUE殘差與各月海域 40.5—41.0 °S、79.0—79.5 °W 位置 對 應(yīng)的 SST、Chl-a、SSTA 和 SSHA 進 行 回歸分析。

根據(jù)式（6）可得基于 SST、Chl-a、SSTA 和SSHA的剩余產(chǎn)量模型EDSP5的最終表達式。

2.3.6 智利竹筴魚EDSP 模型預(yù)測結(jié)果分析 由模型EDSP1、EDSP2、EDSP3、EDSP4、EDSP5 得出預(yù)測的智利竹筴魚漁獲量與實際漁獲量（圖3），可以看出，基于環(huán)境因子 SST、Chl-a、SSTA 和SSHA 的智利竹筴魚EDSP5 模型對于漁獲量的預(yù)測值與實際值最為接近。

圖3 2012—2018年4—8月智利竹筴魚漁業(yè)實際漁獲量與不同模型預(yù)測漁獲量Fig.3 Actual and predicted catches of Chilean jack mackerel by different model from April to August，2012 to 2018

2.3.7 EDSP 模型擬合結(jié)果評價 對5個模型所預(yù)測的2012—2018年4—8月漁獲量結(jié)果使用AIC 值進行檢驗，AIC 值越小表明擬合效果越優(yōu)，MAE、MRE 越小表明對漁獲量的預(yù)測越準確。從表2 可以看出，這5個模型AIC 值均較小，表明5個模型均可對漁獲量進行較為精確的預(yù)測，但剩余產(chǎn)量模型 EDSP5 的 AIC、MAE、MRE 值均最小，分別為51.79、0.86、25.60，表明該模型的擬合度最好且預(yù)測精度最高，為最優(yōu)EDSP 模型，其次為EDSP4，EDSP1 的擬合度相對較差。EDSP4、EDSP5 模型結(jié)果擬合程度相差結(jié)果很小，且對漁獲量的預(yù)測情況差異不大，考慮到建立模型時的方便性與簡潔性，可以選擇環(huán)境因子SST、Chl-a、SSTA來建立智利竹筴魚資源預(yù)測模型。

表2 模型擬合結(jié)果Table 2 Summary results of model

3 討論

本研究利用GAM 模型對東南太平洋智利竹筴魚CPUE 進行標準化時，不僅考慮空間因素對智利竹筴魚CPUE 的影響，也將環(huán)境因子加入解釋變量中，去除空間和環(huán)境雙重因素對智利竹筴魚CPUE的影響，以便更好地構(gòu)建EDSP模型。

從預(yù)測結(jié)果最優(yōu)剩余產(chǎn)量模型EDSP5 模型的表達式可看出，決定智利竹筴魚漁獲量的主要因素除 CPUE 量外，SST、Chl-a、SSTA、SSHA 也應(yīng)考慮其中，SST 是導(dǎo)致智利竹筴魚資源變動的重要海洋環(huán)境因子[3-4]。當(dāng)對應(yīng)月份海域40.5—41.0 °S、79.0—79.5 °W位置的SST低于14 ℃時，智利竹筴魚CPUE 殘差小于0；當(dāng)對應(yīng)月份SST 在14～16 ℃，智利竹筴魚CPUE 殘差大于0，由此可見，當(dāng)海面溫度在14～16 ℃時，海水溫度適宜，該環(huán)境更適合智利竹筴魚的生長，智利竹筴魚資源量會有所增加，在相同捕撈努力量情況下，漁獲量會增加；相反，環(huán)境溫度不適宜時資源量下降，漁獲量減少。該結(jié)果也與牛明香等[4]的研究結(jié)果一致。SST 對智利竹筴魚CPUE 的影響可直接反映到其資源豐度上，進而導(dǎo)致智利竹筴魚漁獲量的改變。研究結(jié)果表明用對應(yīng)月份40.5—41.0 °S、79.0—79.5 °W 海域環(huán)境因子作為表征智利竹筴魚漁場的環(huán)境情況是可行的，且使用海表面環(huán)境因子數(shù)據(jù)可以對智利竹筴魚漁獲量進行研究。

DSP 模型在加入海面溫度數(shù)據(jù)后，假定捕撈努力量值為常數(shù)，智利竹筴魚CPUE與SST存在線性關(guān)系；在加入環(huán)境因子SST 和Chl-a 影響后，發(fā)現(xiàn)2 種環(huán)境因子對智利竹筴魚CPUE 的影響比只考慮海面溫度時的影響大，對智利竹筴魚漁獲量的預(yù)測結(jié)果也比EDSP2 模型中的平均絕對誤差和相對誤差小；EDSP5 模型將4個環(huán)境因子均加入其中，結(jié)果表明其為預(yù)測效果最佳的模型。

智利竹筴魚為中上層海洋魚類，資源容易隨著環(huán)境和氣候變化而發(fā)生變化[28]。本研究表明，SST、Chl-a、SSTA、SSHA對智利竹筴魚資源產(chǎn)生了一定的影響，導(dǎo)致CPUE 的波動。SST是導(dǎo)致智利竹筴魚CPUE 產(chǎn)生月間波動的最重要環(huán)境因子，其次為SSTA、Chl-a，SSHA 則影響相對較小，且隨著有效環(huán)境因子加入得越多，對于漁獲量的分析和預(yù)測就愈加準確和全面，但從EDSP4 與EDSP5模型的結(jié)果對比來看，擬合程度相似，平均相對誤差和平均絕對誤差相差較小，對漁獲量的預(yù)測情況差異不大，考慮到模型的簡潔性及有效性，建立模型時可舍棄海面高度異常值這一環(huán)境因子。

21 世紀初，由于環(huán)境變化和人為捕撈導(dǎo)致資源驟減，為了加強漁業(yè)資源的有效管理，澳大利亞、中國、智利、俄羅斯等國家于2006年共同發(fā)起籌建南太平洋區(qū)域漁業(yè)管理組織計劃，經(jīng)過多方面的交流溝通在2012年實施《南太平洋公海漁業(yè)資源養(yǎng)護和管理公約》，限制各國的捕撈水平，旨在恢復(fù)智利竹筴魚資源量。在時間年際變化趨勢上，2014—2018年智利竹筴魚資源量相比較2012—2013年資源量有所增加，說明南太平洋區(qū)域漁業(yè)管理組織總捕撈許可量管理制度初見成效。近幾年，東南太平洋智利竹筴魚資源一直處于恢復(fù)時期，因此對于智利竹筴魚漁場及資源的研究就尤為重要，本文為未來漁業(yè)管理的開發(fā)利用提供參考依據(jù)，如采取更為科學(xué)嚴格的管理措施，制定相應(yīng)合理的捕撈力度。

剩余產(chǎn)量模型由于對數(shù)據(jù)要求不高且簡單實用等優(yōu)勢在漁業(yè)中應(yīng)用廣泛，近幾年來智利竹筴魚受到環(huán)境變化和人為捕撈影響而產(chǎn)量較低，智利竹筴魚的數(shù)據(jù)量可以滿足該模型的要求。相比較于傳統(tǒng)的剩余產(chǎn)量模型，本文選取了SST、Chl-a、SSTA、SSHA 4個環(huán)境因子，將這些環(huán)境因子加入了所建立的剩余產(chǎn)量模型中，動態(tài)研究和預(yù)測基于環(huán)境變化下智利竹筴魚漁獲量的變化。與李綱等[15]、鄭曉瓊等[16]、易煒等[29]、楊嘉樑等[30]構(gòu)建的剩余產(chǎn)量模型相比，本文考慮了多個海洋環(huán)境因子對資源量的影響，建立基于不同海洋環(huán)境因子的智利竹筴魚的剩余產(chǎn)量模型，更準確地對智利竹筴魚漁獲量進行預(yù)測。

本文只考慮了空間因素和部分環(huán)境因素與智利竹筴魚漁獲量的關(guān)系，在一定程度上會影響預(yù)測結(jié)果；CPUE 只能間接地反映智利竹筴魚資源豐度的變化，本文直接采用CPUE 代表資源豐度可能會使結(jié)果出現(xiàn)誤差；本文僅采用4個環(huán)境因子建立模型，但是對于智利竹筴魚資源量的影響是多方面的，今后的研究也可以加入更多影響漁場資源變動的因子例如溶解氧、海水鹽度等，以便更好、更全面地進行智利竹筴魚資源豐度的研究，為今后的智利竹筴魚資源的開發(fā)提供一定的理論參考。