印度洋北部鸢乌贼CPUE标准化初步研究

温利红,张 衡,方 舟,3,4,5,6*,陈新军,3,4,5,6

(1.上海海洋大学海洋科学学院，上海 201306；2.中国水产科学院东海水产研究所，上海 200090;3.大洋渔业资源可持续开发教育部重点实验室，上海 201306；4.国家远洋渔业工程技术研究中心，上海 201306；5.农业农村部大洋渔业开发重点实验室，上海 201306;6.农业农村部大洋渔业资源环境科学观测实验站，上海 201306)

引 言

鸢乌贼(Sthenoteuthisoualaniensis)，隶属头足纲(Cephalopoda)，柔鱼科(Ommastrephidae)，鸢乌贼属(Sthenoteuthis)。该种类生命周期短，繁殖能力强，资源量丰富，是世界重要的经济性头足类，广泛分布于印度洋、太平洋的赤道和亚热带海域[1-2]。目前印度洋西北部海域和南海是我国捕捞鸢乌贼的主要渔场[3-4]。印度洋沿岸国渔业资源丰富，但除少数发达国家外，大多数国家的经济发展水平不高，渔业资源也尚未得到开发，相对于其他公海鸢乌贼渔场，开发潜力较大，日益受到人们重视[4-5]。由于鸢乌贼具有昼夜垂直迁移习性和趋光性，目前我国针对该物种的主要捕捞方式为灯光敷网、灯光罩网和鱿钓[6-8]。

合理分析环境因素与渔获量的关系，是了解和掌握鸢乌贼资源量及分布规律，合理开发和利用鸢乌贼渔业资源的重要前提之一，但两者之间的关系往往不是简单的线性关系。单位捕捞努力量渔获量(catch per unit effort,CPUE)通常被假定与渔业资源量之间成正比关系，是衡量渔业资源密度的一个重要指标，被广泛应用于渔业资源评估与管理中[9-10]。由于在实际渔业生产过程中，受渔场时空分布、海洋环境、气候变化及作业渔船的捕捞方式等多种因素的影响，CPUE与资源量之间存在较为复杂的关系，不能准确表达出不同作业方式之间的捕捞努力量差异，无法全面、真实地反映渔业资源量及分布情况[11]。在此情况下，需要对CPUE进行标准化处理，去除影响因素，减小误差，从而能够科学的衡量渔业资源量[10]。目前，研究渔业资源与环境之间的关系模型方法有：广义线性模型(generalized linear model,GLM)和广义可加性模型(generalized additive model,GAM)，这两种模型方法已被国内外广泛到渔业中[12]。其中GLM模型主要解决响应变量和解释变量之间的线性问题，由于CPUE受多个变量的影响，且与各影响变量之间的关系较为复杂，可能存在非线性关系，而GAM模型是GLM模型的非参数化拓展，是一种非线性关系模型，能够较好地处理非线性问题，更加有效地探讨渔业资源与影响变量之间的关系，提高研究的准确性[13-14]。为此，本文根据2017—2019年中国远洋渔业协会鱿钓技术组和公海围拖网技术组提供的生产统计数据，利用GLM模型筛选关键时空环境因子，以GAM模型分析印度洋北部鸢乌贼渔场与时空环境因子间的关系，并以此对进行CPUE标准化，为我国对印度洋鸢乌贼渔业资源的合理利用和渔情预报提供相关的参考依据。

1 材料与方法

1.1 材料来源

印度洋北部鸢乌贼生产统计资料来自中国远洋渔业协会鱿钓技术组和公海围拖网技术组，该统计资料包含灯光敷网、灯光罩网和鱿钓3种作业方式。统计内容包括作业日期、作业次数、作业经度、作业纬度和渔获量。空间分辨率为1°×1°，本文研究的海域范围为10°～25°N，50°～75°E(图1)。时间为2017—2019年。

图1 印度洋北部海域范围(注：黑色线框内为研究区域)

海表面温度(Sea Surface Temperature,SST)、海表面高度(Sea Surface Height,SSH)、海表盐度(Sea Surface Salinity,SSS)、光合有效辐射(Photosynthetically Active Radiation,PAR)、风速(Wind Speed,WS)以及流速(Current Speed,U)被认为是影响大洋性鱿鱼渔场分布的重要环境因子[15-16]，因此在本研究中纳入CPUE标准化的影响因子。海表温度和光合有效辐射数据来源于美国大气与海洋局(National Ocean and Atmosphere Administration,NOAA)中太平洋观测网点(https://oceanwatch.pifsc.noaa.gov/erddap/index.html)；海表盐度、海表面高度、风速和流速数据来源于夏威夷大学网站(http://apdrc.soest.hawaili.edu/data/data.php)。时间分辨率为月，空间分辨率为0.5°×0.5°。

1.2 数据预处理

1.2.1 计算

单位捕捞努力量渔获量(catch per unit effort, CPUE)可以作为表征鸢乌贼资源密度的指标，计算公式如下[17]：

(1)

式中，CPUE单位为t/次；C表示一艘渔船一天的产量；E表示其对应的作业次数，灯光敷网和灯光罩网以网次计算、鱿钓是按照每天作业位置的变化次数来计算。

1.2.2 数据匹配

通过Excel和Matlab软件的相关程序包，将下载的海表温度(SST)、海表盐度(SSS)、海表面高度(SSH)、光合有效辐射(PAR)、风速(WS)和流速(U)等环境因子，利用克里金插值法与渔业数据(经纬度、产量、标准化CPUE以及捕捞方式)进行匹配[6-7]，使得环境数据与渔业数据一一对应，均调整为0.5°×0.5°的分辨率。

1.3 变量的相互独立性检验

本研究根据鸢乌贼渔业特点选择解释变量为年、月、经度、纬度、海表面温度、海表盐度、海表面高度、光合有效辐射、风速、流速和作业方式(fishing types, FT)。利用方差膨胀系数(VIF)和Spearman相关系数对解释变量进行相互独立性检验。当VIF<10时，表明变量之间不存在多重共线性，VIF>10时，表明变量之间存在严重的多重共线性[18]。根据独立性检验，结果显示，各个解释变量的VIF均小于10，表明解释变量之间不存在严重多重共线性问题(表1)。

表1 解释变量间方差膨胀因子

1.4 研究方法

目前CPUE标准化方法主要采用广义线性模型(GLM)和广义加性模型(GAM)。国内外较多研究表明，利用GLM模型可获得影响CPUE的主要因子及其贡献度,并可考虑多因子的交互效应；而GAM模型在CPUE标准化中的效果优于GLM模型[12,19-20]。因此本研究首先利用广义线性模型(GLM)对影响鸢乌贼资源密度的因子进行显著性检验，筛选出具有显著性影响的因子，随后利用广义加性模型(GAM)对显著性因子和鸢乌贼资源密度关系进行研究，并利用该模型尝试对不同时间尺度(年和月)鸢乌贼CPUE进行标准化。

GLM模型假设响应变量的期望值与解释变量呈线性关系[21]：

(2)

式中，g为链接函数;μi=E(Yi),Yi为第i个响应变量;Xi为第i个响应变量的解释变量;β为模型估计参数。本研究假设CPUE服从对数正态分布，则GLM模型表达式为：

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

(3)

式中，CPUE为每艘船每次作业的捕捞产量，对CPUE加上常数1是为解决CPUE出现零值的现象；ε为u误差项，假设其服从正态分布。GLM模型中，将时间(年、月)、空间(经度、纬度)、环境(SST、SSS、PAR、WS、SSH、U)和作业方式因子作为解释变量，其中变量年、月、经度、纬度和作业方式作为离散变量，其他变量(SST、SSS、PAR、WS、SSH、U)作为连续变量。

GAM是GLM模型的非线性拓展，通过将函数与相应变量进行移动的变化，将基于指数分布的回归与一般线性回归进行整合[22]，即：

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

(4)

式中，g为链接函数;μi=E(Yi)，Yi为第i个响应变量;Xi为第i个响应变量的解释变量;ε为模型估计参数;fi为平滑函数。根据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是为防止响应变量出现零值，再进行了对数化处理；s为自然立方样条平滑(natural cube spline smoother)；s(year)为年效应；s(month)为月效应；s(latitude)为纬度效应；s(longitude)为经度效应；s(SST)为海表温度效应；s(SSS)为海表盐度效应；s(PAR)为光合有效辐射效应；s(WS)为风速效应；s(U)为流速效应；由于作业方式的数值组成太少，因此将作业方式以因子factor(FT)形式处理。

根据赤池信息量准则(AIC)值，选取最佳模型[23]。AIC值计算如下：

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

(6)

式中，m为模型中参数的个数。

本研究使用R(V3.2.2)进行统计分析与处理。

2 结果与分析

2.1 解释变量ln(CPUE+1)的统计分布检验

首先验证ln(CPUE+1)是否服从正态分布，经K-S检验，ln(CPUE+1)的数据点在正态Q-Q图中基本形成一条直线(图2-B)，ln(CPUE+1)服从正态分布(图2-A)，这说明本文研究中关于ln(CPUE+1)服从正态分布的假设是合理的，可以运用GLM模型和GAM模型进行分析。

图2 2017—2019年印度洋北部鸢乌贼ln(CPUE+1)的频次分布及其检验

2.2 模型分析

2.2.1 GLM模型分析

通过GLM模型对各因子进行显著性检验(见表2)。本文以P<0.05来确认变量是否具有显著性，结果表明，年、月、纬度、SST、WS、U以及作业方式均为显著性变量，且除因子SST对CPUE的影响为显著性外，其他因子都为极显著性(P<0.01)；经度、SSS、PAR及SSH为不显著性变量(P>0.05)，对CPUE的影响不明显。因此，选择7个显著性解释变量(年、月、纬度、SST、WS、U以及作业方式)纳入GAM模型对CPUE进行标准化。

表2 GLM模型自变量显著性检验

2.2.2 GAM模型分析

随后将上述具有显著性的解释变量年、月、纬度、SST、WS、U以及作业方式逐一加入GAM模型中，进行运算分析，发现所选取的因子建立的GAM模型，AIC值最小，其拟合效果最好(见表3)，最终得到最佳的GAM模型为：

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

(7)

并且根据P值，得出年、月、纬度、SST、WS、U及作业方式均为显著性变量，对CPUE的影响都为极显著性(P<0.01)。模型对CPUE的总偏差解释为21.6%，其中作业方式变量对CPUE的影响最大，解释了5.6%的总偏差，说明作业方式对CPUE的影响最大；随后影响由大到小依次是年(5.22%)、纬度(4.94%)、月(3.94%)、SST(1.2%)、U(0.6%)、WS(0.1%)。

表3 GAM模型统计结果

2.3 各因子对CPUE的影响

从时间因子来看，年为极显著性变量，对CPUE具有较大影响；CPUE总体在2017—2018年呈增长趋势，至2018年达到峰值，随后呈小幅度下降趋势，2019年较为集中。月变化对CPUE的影响也极显著，总体变化幅度较大，具有明显的季节变化，其中1—7月呈下降趋势，并在7月达到最低值，此后随月份逐渐上升，并在12月份达到最大值。从空间因子来看，纬度为极显著性变量，变化趋势总体呈波浪式起伏，且变化幅度较大，首先在10°N～12°N处呈下降趋势，并在12°N附近达最小值，后随着纬度的增加呈上升趋势，于18°N附近达到最大值，随后又呈下降趋势，其中CPUE在14°N～17°N范围分布较为集中(图3)。

图3 基于GAM模型的时空与环境因子效应对印度洋北部鸢乌贼CPUE的影响

从环境因子来看，SST为极显著性变量，由图2可知，印度洋北部鸢乌贼作业区域的SST分布范围为23～31℃，CPUE主要集中分布在26～28.5 ℃之间，CPUE在23～28 ℃范围内呈小幅度下降趋势，至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；其余时间段内标准化CPUE均高于或接近对应的名义CPUE。其中2018年11月，名义CPUE和标准化CPUE都达到最大值；2017年2月，名义CPUE达到最小值，2018年7月，标准化CPUE达最小值。总体上名义CPUE的变化波动幅度较大，而标准化CPUE变化幅度较小，且名义CPUE与标准化后的CPUE变化趋势基本一致(图4b)。

图4 基于GAM模型的印度洋北部鸢乌贼年(a)和月(b)平均标准化CPUE

3 讨论与分析

3.1 时空因子对CPUE的影响

从时间因子来看，GAM模型结果表明，年对CPUE的贡献率为5.22%。其中，2018年CPUE的影响效应比其余年份高，主要是因为2018年的生产总量和频次数据相对较多，因此所统计的总产量较其他年份多，对资源密影响相对较大。同时月间波动性较大，10月到次年1月对CPUE的影响较大。鸢乌贼为一年生群体，全年可捕捞作业，一方面，夏季受西南季风的高强度且持久性的影响，索马里沿岸上升流不断将海底丰富的营养盐带到海表层，使得浮游动植物大量繁殖，当上升流到达索马里半岛外的大涡旋处，在季风和大涡旋的共同作用下产生强大的离岸流，将近海岸丰富的营养盐和叶绿素以及沿海独特的浮游植物运输到阿拉伯海中部海域，此过程存在一定的滞后性，所以秋冬季季节生产力也较高[24-25]；另一方面，9月以后，随着夏季风的消退、相关的海洋和气候环境的变化，包括云量减少、风速减弱、降雨减少等因素的综合影响下，秋冬季的初级生产力和次级生产力显著增加，生物量也显著升高[26]。此外，夏季受西南季风的影响，印度洋北部海域风浪较大，不利于捕捞作业，同时由于2017年6—7月渔场转移其其他海域，导致印度洋北部海域无生产数据，造成该两月对CPUE的影响程度较低。

从空间因子来看，GAM模型结果表明，纬度对CPUE的贡献率为4.94%，空间范围在14°N～18°N、59°E～63°E尤其对CPUE影响较大。该海域受索马里海流和季风海流的影响，营养盐丰富，表层藻类繁盛，叶绿素浓度增高，提供了丰富的饵料，是鸢乌贼的一个重要渔场所在地[27-28]。根据本文统计数据可知，渔船作业位置主要集中在16°N～18°N，61.5°E～63°E海域，该海域是鸢乌贼的中心渔场，资源量丰富[29]。陈新军等[30]研究也认为印度洋鸢乌贼渔场在12°N～18°N、58°E～61°E具有较高的CPUE，本文研究结果与其结论较一致。

3.2 环境因子对CPUE的影响

在所研究的环境因子中，SST对CPUE的影响最大。鸢乌贼作为暖水性的大洋性头足类，其渔场分布与SST有着高度的关联性[31]。陈新军等[3]研究认为，印度洋中心渔场的水温为26.4～29.0 ℃；余为等[32]在印度洋西北海域鸢乌贼的栖息地研究中，得出鸢乌贼的最适SST范围为27～29 ℃。本文研究中，作业温度集中在26～28.5 ℃，与上述研究结果较一致。与此同时，较大尺度气候因子也影响着SST的变化[33]。其中印度洋偶极子(Indian Ocean Dipole,IOD)是影响印度洋气候环境变化的重要特殊气候现象：当IOD发生时，印度洋处于海温异常状态，主要通过海—气相互作用影响印度洋海盆海洋环境，IOD处于正相位时期，西印度洋海温高于东印度洋，且印度洋季风环流减弱；相反，处于负相位时期，西印度洋海温低于东印度洋，季风环流增强[34]。Wei等[35]研究表明，赤道印度洋海域的印度洋偶极子指数(DMI)从2017—2019年经历了正负波动到最大正值的过程，IOD从负相位转变为正相位，西印度洋海表面温度逐渐升高。同时Anderson等[36]的研究结果也表明这种偶极子分布形态存在于海表温度的变化中。

另外海流和风速对印度洋鸢乌贼CPUE也有着较大影响，这在本研究的GAM模型结果里也有所体现。印度洋北部被非洲大陆、亚洲大陆三面环抱，使得海洋环境和气候变化深受大陆的影响，成为世界典型的季风海域，季风是驱动印度洋环境变化和海洋生态动力过程的最主要因素之一[37]。印度洋北部受气压带和风带位置的季节性移动，以及海陆热力性质差异的影响，海域多为风成海流和上升流，进而海域表层环流促使深层和表层的海水进行交换，将海水深层丰富的营养盐带到海表层，带来较高的初级生产力，使得浮游动植物生物量大量富集，并呈现出明显的区域差异和季节变化[38]，这些因素都会对鸢乌贼资源量和空间分布产生潜在影响。此外，由于作业方式的差异，风速的大小对鱿钓作业方式影响很大，风速大会使钓具上下浮动频率变大，影响捕捞效率[39]。但对灯光敷网和灯光罩网捕捞作业影响较小。全球气候变化对印度洋海洋环境变化产生重要影响，影响的因素也是多方面的。因此，可以考虑将U和WS作为突出的潜在因子应用于后续的渔情预报中。

3.3 不同作业方式对CPUE的影响

以往的研究中，由于作业方式较为单一，往往会忽略其他因素对CPUE所产生的影响，而在印度洋鸢乌贼渔业中，存在三种作业方式，不同作业方式的CPUE与单位渔船数量、单位渔船主机功率、渔具作业性能、捕捞作业技术水平和作业时间等因素有关[39-40]。本研究中灯光敷网和灯光罩网作业方式较相似，因此两者之间差异较小。然而，鱿钓作业对CPUE的影响与前两者相比有着较大的差异，GAM模型中的结果也说明了这一问题。灯光敷网和灯光罩网为主动捕捞，主要的捕捞对象是趋光性的鱿鱼和中上层鱼类，渔获选择性较低，作业位置转移较为灵活，捕捞努力量较高，且捕捞方式受到海洋环境变化的影响较小，因此捕捞效率的变化较小，CPUE起伏波动较小；而鱿钓的作业方式相对被动，且捕捞位置相对固定，渔获选择性高，易受捕捞对象资源量变动的制约，且钓具易受海流、风浪的影响，因此导致捕捞效率变化较大，CPUE起伏波动较大。因此，后续的相关研究中，不仅需要考虑时空环境因子对印度洋鸢乌贼CPUE的影响，也不应忽视捕捞方式所带来的影响。

4 结论

本文根据2017—2019年印度洋北部鸢乌贼渔业生产数据，通过使用GLM和GAM模型，分析时空因子和重要的环境因子对印度洋北部鸢乌贼CPUE变化的影响，以期为科学管理和利用鸢乌贼渔业资源提供理论依据。在后续的研究中，应考虑因子间的综合效应以及加强时间序列的样本采集工作；同时在比较不同作业方式的差异时，还应该考虑到作业方式对资源量反映的准确性，结合多种方法综合考虑资源量和环境因子间的关系。