印度洋北部鳶烏賊CPUE標準化初步研究

溫利紅,張 衡,方 舟,3,4,5,6*,陳新軍,3,4,5,6

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

引 言

鳶烏賊(Sthenoteuthisoualaniensis)，隸屬頭足綱(Cephalopoda)，柔魚科(Ommastrephidae)，鳶烏賊屬(Sthenoteuthis)。該種類生命周期短，繁殖能力強，資源量豐富，是世界重要的經(jīng)濟性頭足類，廣泛分布于印度洋、太平洋的赤道和亞熱帶海域[1-2]。目前印度洋西北部海域和南海是我國捕撈鳶烏賊的主要漁場[3-4]。印度洋沿岸國漁業(yè)資源豐富，但除少數(shù)發(fā)達國家外，大多數(shù)國家的經(jīng)濟發(fā)展水平不高，漁業(yè)資源也尚未得到開發(fā)，相對于其他公海鳶烏賊漁場，開發(fā)潛力較大，日益受到人們重視[4-5]。由于鳶烏賊具有晝夜垂直遷移習性和趨光性，目前我國針對該物種的主要捕撈方式為燈光敷網(wǎng)、燈光罩網(wǎng)和魷釣[6-8]。

合理分析環(huán)境因素與漁獲量的關系，是了解和掌握鳶烏賊資源量及分布規(guī)律，合理開發(fā)和利用鳶烏賊漁業(yè)資源的重要前提之一，但兩者之間的關系往往不是簡單的線性關系。單位捕撈努力量漁獲量(catch per unit effort,CPUE)通常被假定與漁業(yè)資源量之間成正比關系，是衡量漁業(yè)資源密度的一個重要指標，被廣泛應用于漁業(yè)資源評估與管理中[9-10]。由于在實際漁業(yè)生產(chǎn)過程中，受漁場時空分布、海洋環(huán)境、氣候變化及作業(yè)漁船的捕撈方式等多種因素的影響，CPUE與資源量之間存在較為復雜的關系，不能準確表達出不同作業(yè)方式之間的捕撈努力量差異，無法全面、真實地反映漁業(yè)資源量及分布情況[11]。在此情況下，需要對CPUE進行標準化處理，去除影響因素，減小誤差，從而能夠科學的衡量漁業(yè)資源量[10]。目前，研究漁業(yè)資源與環(huán)境之間的關系模型方法有：廣義線性模型(generalized linear model,GLM)和廣義可加性模型(generalized additive model,GAM)，這兩種模型方法已被國內(nèi)外廣泛到漁業(yè)中[12]。其中GLM模型主要解決響應變量和解釋變量之間的線性問題，由于CPUE受多個變量的影響，且與各影響變量之間的關系較為復雜，可能存在非線性關系，而GAM模型是GLM模型的非參數(shù)化拓展，是一種非線性關系模型，能夠較好地處理非線性問題，更加有效地探討漁業(yè)資源與影響變量之間的關系，提高研究的準確性[13-14]。為此，本文根據(jù)2017—2019年中國遠洋漁業(yè)協(xié)會魷釣技術(shù)組和公海圍拖網(wǎng)技術(shù)組提供的生產(chǎn)統(tǒng)計數(shù)據(jù)，利用GLM模型篩選關鍵時空環(huán)境因子，以GAM模型分析印度洋北部鳶烏賊漁場與時空環(huán)境因子間的關系，并以此對進行CPUE標準化，為我國對印度洋鳶烏賊漁業(yè)資源的合理利用和漁情預報提供相關的參考依據(jù)。

1 材料與方法

1.1 材料來源

印度洋北部鳶烏賊生產(chǎn)統(tǒng)計資料來自中國遠洋漁業(yè)協(xié)會魷釣技術(shù)組和公海圍拖網(wǎng)技術(shù)組，該統(tǒng)計資料包含燈光敷網(wǎng)、燈光罩網(wǎng)和魷釣3種作業(yè)方式。統(tǒng)計內(nèi)容包括作業(yè)日期、作業(yè)次數(shù)、作業(yè)經(jīng)度、作業(yè)緯度和漁獲量。空間分辨率為1°×1°，本文研究的海域范圍為10°～25°N，50°～75°E(圖1)。時間為2017—2019年。

圖1 印度洋北部海域范圍(注：黑色線框內(nèi)為研究區(qū)域)

海表面溫度(Sea Surface Temperature,SST)、海表面高度(Sea Surface Height,SSH)、海表鹽度(Sea Surface Salinity,SSS)、光合有效輻射(Photosynthetically Active Radiation,PAR)、風速(Wind Speed,WS)以及流速(Current Speed,U)被認為是影響大洋性魷魚漁場分布的重要環(huán)境因子[15-16]，因此在本研究中納入CPUE標準化的影響因子。海表溫度和光合有效輻射數(shù)據(jù)來源于美國大氣與海洋局(National Ocean and Atmosphere Administration,NOAA)中太平洋觀測網(wǎng)點(https://oceanwatch.pifsc.noaa.gov/erddap/index.html)；海表鹽度、海表面高度、風速和流速數(shù)據(jù)來源于夏威夷大學網(wǎng)站(http://apdrc.soest.hawaili.edu/data/data.php)。時間分辨率為月，空間分辨率為0.5°×0.5°。

1.2 數(shù)據(jù)預處理

1.2.1 計算

單位捕撈努力量漁獲量(catch per unit effort, CPUE)可以作為表征鳶烏賊資源密度的指標，計算公式如下[17]：

(1)

式中，CPUE單位為t/次；C表示一艘漁船一天的產(chǎn)量；E表示其對應的作業(yè)次數(shù)，燈光敷網(wǎng)和燈光罩網(wǎng)以網(wǎng)次計算、魷釣是按照每天作業(yè)位置的變化次數(shù)來計算。

1.2.2 數(shù)據(jù)匹配

通過Excel和Matlab軟件的相關程序包，將下載的海表溫度(SST)、海表鹽度(SSS)、海表面高度(SSH)、光合有效輻射(PAR)、風速(WS)和流速(U)等環(huán)境因子，利用克里金插值法與漁業(yè)數(shù)據(jù)(經(jīng)緯度、產(chǎn)量、標準化CPUE以及捕撈方式)進行匹配[6-7]，使得環(huán)境數(shù)據(jù)與漁業(yè)數(shù)據(jù)一一對應，均調(diào)整為0.5°×0.5°的分辨率。

1.3 變量的相互獨立性檢驗

本研究根據(jù)鳶烏賊漁業(yè)特點選擇解釋變量為年、月、經(jīng)度、緯度、海表面溫度、海表鹽度、海表面高度、光合有效輻射、風速、流速和作業(yè)方式(fishing types, FT)。利用方差膨脹系數(shù)(VIF)和Spearman相關系數(shù)對解釋變量進行相互獨立性檢驗。當VIF<10時，表明變量之間不存在多重共線性，VIF>10時，表明變量之間存在嚴重的多重共線性[18]。根據(jù)獨立性檢驗，結(jié)果顯示，各個解釋變量的VIF均小于10，表明解釋變量之間不存在嚴重多重共線性問題(表1)。

表1 解釋變量間方差膨脹因子

1.4 研究方法

目前CPUE標準化方法主要采用廣義線性模型(GLM)和廣義加性模型(GAM)。國內(nèi)外較多研究表明，利用GLM模型可獲得影響CPUE的主要因子及其貢獻度,并可考慮多因子的交互效應；而GAM模型在CPUE標準化中的效果優(yōu)于GLM模型[12,19-20]。因此本研究首先利用廣義線性模型(GLM)對影響鳶烏賊資源密度的因子進行顯著性檢驗，篩選出具有顯著性影響的因子，隨后利用廣義加性模型(GAM)對顯著性因子和鳶烏賊資源密度關系進行研究，并利用該模型嘗試對不同時間尺度(年和月)鳶烏賊CPUE進行標準化。

GLM模型假設響應變量的期望值與解釋變量呈線性關系[21]：

(2)

式中，g為鏈接函數(shù);μi=E(Yi),Yi為第i個響應變量;Xi為第i個響應變量的解釋變量;β為模型估計參數(shù)。本研究假設CPUE服從對數(shù)正態(tài)分布，則GLM模型表達式為：

ln(CPUE+1)=Year+Month+Latitude+Longitude+SST+SSS+PAR+WS+U+FT+ε

(3)

式中，CPUE為每艘船每次作業(yè)的捕撈產(chǎn)量，對CPUE加上常數(shù)1是為解決CPUE出現(xiàn)零值的現(xiàn)象；ε為u誤差項，假設其服從正態(tài)分布。GLM模型中，將時間(年、月)、空間(經(jīng)度、緯度)、環(huán)境(SST、SSS、PAR、WS、SSH、U)和作業(yè)方式因子作為解釋變量，其中變量年、月、經(jīng)度、緯度和作業(yè)方式作為離散變量，其他變量(SST、SSS、PAR、WS、SSH、U)作為連續(xù)變量。

GAM是GLM模型的非線性拓展，通過將函數(shù)與相應變量進行移動的變化，將基于指數(shù)分布的回歸與一般線性回歸進行整合[22]，即：

g(μi)=α+∑i=1fi(xi)+ε

(4)

式中，g為鏈接函數(shù);μi=E(Yi)，Yi為第i個響應變量;Xi為第i個響應變量的解釋變量;ε為模型估計參數(shù);fi為平滑函數(shù)。根據(jù)GLM模型篩選出的顯著性解釋變量依次加入GAM模型，則GAM模型的表達式為：

ln(CPUE+1)～s(Year)+s(Month)+s(Latitude)+s(Longitude)+
s(SST)+s(SSS)+s(PAR)+s(WS)+s(U)+factor(FT)+ε

(5)

式中，對CPUE+1是為防止響應變量出現(xiàn)零值，再進行了對數(shù)化處理；s為自然立方樣條平滑(natural cube spline smoother)；s(year)為年效應；s(month)為月效應；s(latitude)為緯度效應；s(longitude)為經(jīng)度效應；s(SST)為海表溫度效應；s(SSS)為海表鹽度效應；s(PAR)為光合有效輻射效應；s(WS)為風速效應；s(U)為流速效應；由于作業(yè)方式的數(shù)值組成太少，因此將作業(yè)方式以因子factor(FT)形式處理。

根據(jù)赤池信息量準則(AIC)值，選取最佳模型[23]。AIC值計算如下：

AIC=-2lnl(p1,p2,…,pm,σ2)+2m

(6)

式中，m為模型中參數(shù)的個數(shù)。

本研究使用R(V3.2.2)進行統(tǒng)計分析與處理。

2 結(jié)果與分析

2.1 解釋變量ln(CPUE+1)的統(tǒng)計分布檢驗

首先驗證ln(CPUE+1)是否服從正態(tài)分布，經(jīng)K-S檢驗，ln(CPUE+1)的數(shù)據(jù)點在正態(tài)Q-Q圖中基本形成一條直線(圖2-B)，ln(CPUE+1)服從正態(tài)分布(圖2-A)，這說明本文研究中關于ln(CPUE+1)服從正態(tài)分布的假設是合理的，可以運用GLM模型和GAM模型進行分析。

圖2 2017—2019年印度洋北部鳶烏賊ln(CPUE+1)的頻次分布及其檢驗

2.2 模型分析

2.2.1 GLM模型分析

通過GLM模型對各因子進行顯著性檢驗(見表2)。本文以P<0.05來確認變量是否具有顯著性，結(jié)果表明，年、月、緯度、SST、WS、U以及作業(yè)方式均為顯著性變量，且除因子SST對CPUE的影響為顯著性外，其他因子都為極顯著性(P<0.01)；經(jīng)度、SSS、PAR及SSH為不顯著性變量(P>0.05)，對CPUE的影響不明顯。因此，選擇7個顯著性解釋變量(年、月、緯度、SST、WS、U以及作業(yè)方式)納入GAM模型對CPUE進行標準化。

表2 GLM模型自變量顯著性檢驗

2.2.2 GAM模型分析

隨后將上述具有顯著性的解釋變量年、月、緯度、SST、WS、U以及作業(yè)方式逐一加入GAM模型中，進行運算分析，發(fā)現(xiàn)所選取的因子建立的GAM模型，AIC值最小，其擬合效果最好(見表3)，最終得到最佳的GAM模型為：

ln(CPUE+1)～s(Year)+s(Month)+s(Latitude)+s(SST)+s(WS)+s(U)+factor(FT)

(7)

并且根據(jù)P值，得出年、月、緯度、SST、WS、U及作業(yè)方式均為顯著性變量，對CPUE的影響都為極顯著性(P<0.01)。模型對CPUE的總偏差解釋為21.6%，其中作業(yè)方式變量對CPUE的影響最大，解釋了5.6%的總偏差，說明作業(yè)方式對CPUE的影響最大；隨后影響由大到小依次是年(5.22%)、緯度(4.94%)、月(3.94%)、SST(1.2%)、U(0.6%)、WS(0.1%)。

表3 GAM模型統(tǒng)計結(jié)果

2.3 各因子對CPUE的影響

從時間因子來看，年為極顯著性變量，對CPUE具有較大影響；CPUE總體在2017—2018年呈增長趨勢，至2018年達到峰值，隨后呈小幅度下降趨勢，2019年較為集中。月變化對CPUE的影響也極顯著，總體變化幅度較大，具有明顯的季節(jié)變化，其中1—7月呈下降趨勢，并在7月達到最低值，此后隨月份逐漸上升，并在12月份達到最大值。從空間因子來看，緯度為極顯著性變量，變化趨勢總體呈波浪式起伏，且變化幅度較大，首先在10°N～12°N處呈下降趨勢，并在12°N附近達最小值，后隨著緯度的增加呈上升趨勢，于18°N附近達到最大值，隨后又呈下降趨勢，其中CPUE在14°N～17°N范圍分布較為集中(圖3)。

圖3 基于GAM模型的時空與環(huán)境因子效應對印度洋北部鳶烏賊CPUE的影響

從環(huán)境因子來看，SST為極顯著性變量，由圖2可知，印度洋北部鳶烏賊作業(yè)區(qū)域的SST分布范圍為23～31℃，CPUE主要集中分布在26～28.5 ℃之間，CPUE在23～28 ℃范圍內(nèi)呈小幅度下降趨勢，至28 ℃達到最小，隨后逐漸上升，于30.5 ℃達到最大值。CPUE隨WS增大而緩慢下降趨勢，主要集中分布在5～7 m/s之間。CPUE隨U的增大呈波浪式下降趨勢，在0.05 m/s附近達到最大值，主要集中分布在0.01～0.05 m/s之間(圖3)。

2.4 名義CPUE和標準化CPUE比較

從圖4可知，名義CPUE和基于GAM模型標準化處理后的CPUE變化趨勢隨著年份的增加呈上升趨勢，于2019年達到最高值；整體上，標準化CPUE均略高于名義CPUE，且標準化后的CPUE變化趨勢與名義CPUE變化趨勢一致(圖4a)。

2017年9—12月，2018年4—5月、9月、11—12月，2019年1—4月名義CPUE高于標準化CPUE；其余時間段內(nèi)標準化CPUE均高于或接近對應的名義CPUE。其中2018年11月，名義CPUE和標準化CPUE都達到最大值；2017年2月，名義CPUE達到最小值，2018年7月，標準化CPUE達最小值?？傮w上名義CPUE的變化波動幅度較大，而標準化CPUE變化幅度較小，且名義CPUE與標準化后的CPUE變化趨勢基本一致(圖4b)。

圖4 基于GAM模型的印度洋北部鳶烏賊年(a)和月(b)平均標準化CPUE

3 討論與分析

3.1 時空因子對CPUE的影響

從時間因子來看，GAM模型結(jié)果表明，年對CPUE的貢獻率為5.22%。其中，2018年CPUE的影響效應比其余年份高，主要是因為2018年的生產(chǎn)總量和頻次數(shù)據(jù)相對較多，因此所統(tǒng)計的總產(chǎn)量較其他年份多，對資源密影響相對較大。同時月間波動性較大，10月到次年1月對CPUE的影響較大。鳶烏賊為一年生群體，全年可捕撈作業(yè)，一方面，夏季受西南季風的高強度且持久性的影響，索馬里沿岸上升流不斷將海底豐富的營養(yǎng)鹽帶到海表層，使得浮游動植物大量繁殖，當上升流到達索馬里半島外的大渦旋處，在季風和大渦旋的共同作用下產(chǎn)生強大的離岸流，將近海岸豐富的營養(yǎng)鹽和葉綠素以及沿海獨特的浮游植物運輸?shù)桨⒗Ｖ胁亢Ｓ?，此過程存在一定的滯后性，所以秋冬季季節(jié)生產(chǎn)力也較高[24-25]；另一方面，9月以后，隨著夏季風的消退、相關的海洋和氣候環(huán)境的變化，包括云量減少、風速減弱、降雨減少等因素的綜合影響下，秋冬季的初級生產(chǎn)力和次級生產(chǎn)力顯著增加，生物量也顯著升高[26]。此外，夏季受西南季風的影響，印度洋北部海域風浪較大，不利于捕撈作業(yè)，同時由于2017年6—7月漁場轉(zhuǎn)移其其他海域，導致印度洋北部海域無生產(chǎn)數(shù)據(jù)，造成該兩月對CPUE的影響程度較低。

從空間因子來看，GAM模型結(jié)果表明，緯度對CPUE的貢獻率為4.94%，空間范圍在14°N～18°N、59°E～63°E尤其對CPUE影響較大。該海域受索馬里海流和季風海流的影響，營養(yǎng)鹽豐富，表層藻類繁盛，葉綠素濃度增高，提供了豐富的餌料，是鳶烏賊的一個重要漁場所在地[27-28]。根據(jù)本文統(tǒng)計數(shù)據(jù)可知，漁船作業(yè)位置主要集中在16°N～18°N，61.5°E～63°E海域，該海域是鳶烏賊的中心漁場，資源量豐富[29]。陳新軍等[30]研究也認為印度洋鳶烏賊漁場在12°N～18°N、58°E～61°E具有較高的CPUE，本文研究結(jié)果與其結(jié)論較一致。

3.2 環(huán)境因子對CPUE的影響

在所研究的環(huán)境因子中，SST對CPUE的影響最大。鳶烏賊作為暖水性的大洋性頭足類，其漁場分布與SST有著高度的關聯(lián)性[31]。陳新軍等[3]研究認為，印度洋中心漁場的水溫為26.4～29.0 ℃；余為等[32]在印度洋西北海域鳶烏賊的棲息地研究中，得出鳶烏賊的最適SST范圍為27～29 ℃。本文研究中，作業(yè)溫度集中在26～28.5 ℃，與上述研究結(jié)果較一致。與此同時，較大尺度氣候因子也影響著SST的變化[33]。其中印度洋偶極子(Indian Ocean Dipole,IOD)是影響印度洋氣候環(huán)境變化的重要特殊氣候現(xiàn)象：當IOD發(fā)生時，印度洋處于海溫異常狀態(tài)，主要通過海—氣相互作用影響印度洋海盆海洋環(huán)境，IOD處于正相位時期，西印度洋海溫高于東印度洋，且印度洋季風環(huán)流減弱；相反，處于負相位時期，西印度洋海溫低于東印度洋，季風環(huán)流增強[34]。Wei等[35]研究表明，赤道印度洋海域的印度洋偶極子指數(shù)(DMI)從2017—2019年經(jīng)歷了正負波動到最大正值的過程，IOD從負相位轉(zhuǎn)變?yōu)檎辔?，西印度洋海表面溫度逐漸升高。同時Anderson等[36]的研究結(jié)果也表明這種偶極子分布形態(tài)存在于海表溫度的變化中。

另外海流和風速對印度洋鳶烏賊CPUE也有著較大影響，這在本研究的GAM模型結(jié)果里也有所體現(xiàn)。印度洋北部被非洲大陸、亞洲大陸三面環(huán)抱，使得海洋環(huán)境和氣候變化深受大陸的影響，成為世界典型的季風海域，季風是驅(qū)動印度洋環(huán)境變化和海洋生態(tài)動力過程的最主要因素之一[37]。印度洋北部受氣壓帶和風帶位置的季節(jié)性移動，以及海陸熱力性質(zhì)差異的影響，海域多為風成海流和上升流，進而海域表層環(huán)流促使深層和表層的海水進行交換，將海水深層豐富的營養(yǎng)鹽帶到海表層，帶來較高的初級生產(chǎn)力，使得浮游動植物生物量大量富集，并呈現(xiàn)出明顯的區(qū)域差異和季節(jié)變化[38]，這些因素都會對鳶烏賊資源量和空間分布產(chǎn)生潛在影響。此外，由于作業(yè)方式的差異，風速的大小對魷釣作業(yè)方式影響很大，風速大會使釣具上下浮動頻率變大，影響捕撈效率[39]。但對燈光敷網(wǎng)和燈光罩網(wǎng)捕撈作業(yè)影響較小。全球氣候變化對印度洋海洋環(huán)境變化產(chǎn)生重要影響，影響的因素也是多方面的。因此，可以考慮將U和WS作為突出的潛在因子應用于后續(xù)的漁情預報中。

3.3 不同作業(yè)方式對CPUE的影響

以往的研究中，由于作業(yè)方式較為單一，往往會忽略其他因素對CPUE所產(chǎn)生的影響，而在印度洋鳶烏賊漁業(yè)中，存在三種作業(yè)方式，不同作業(yè)方式的CPUE與單位漁船數(shù)量、單位漁船主機功率、漁具作業(yè)性能、捕撈作業(yè)技術(shù)水平和作業(yè)時間等因素有關[39-40]。本研究中燈光敷網(wǎng)和燈光罩網(wǎng)作業(yè)方式較相似，因此兩者之間差異較小。然而，魷釣作業(yè)對CPUE的影響與前兩者相比有著較大的差異，GAM模型中的結(jié)果也說明了這一問題。燈光敷網(wǎng)和燈光罩網(wǎng)為主動捕撈，主要的捕撈對象是趨光性的魷魚和中上層魚類，漁獲選擇性較低，作業(yè)位置轉(zhuǎn)移較為靈活，捕撈努力量較高，且捕撈方式受到海洋環(huán)境變化的影響較小，因此捕撈效率的變化較小，CPUE起伏波動較??；而魷釣的作業(yè)方式相對被動，且捕撈位置相對固定，漁獲選擇性高，易受捕撈對象資源量變動的制約，且釣具易受海流、風浪的影響，因此導致捕撈效率變化較大，CPUE起伏波動較大。因此，后續(xù)的相關研究中，不僅需要考慮時空環(huán)境因子對印度洋鳶烏賊CPUE的影響，也不應忽視捕撈方式所帶來的影響。

4 結(jié)論

本文根據(jù)2017—2019年印度洋北部鳶烏賊漁業(yè)生產(chǎn)數(shù)據(jù)，通過使用GLM和GAM模型，分析時空因子和重要的環(huán)境因子對印度洋北部鳶烏賊CPUE變化的影響，以期為科學管理和利用鳶烏賊漁業(yè)資源提供理論依據(jù)。在后續(xù)的研究中，應考慮因子間的綜合效應以及加強時間序列的樣本采集工作；同時在比較不同作業(yè)方式的差異時，還應該考慮到作業(yè)方式對資源量反映的準確性，結(jié)合多種方法綜合考慮資源量和環(huán)境因子間的關系。