Study on fishery resource assessment of Pacific saury by length-based cohort analysis
-
摘要:
秋刀鱼 (Cololabis saira) 分布于西北太平洋亚热带到温带海域，是中国远洋渔业主要的捕捞对象之一。为探究其资源状况，根据2014—2018年西北太平洋秋刀鱼的渔获体长组成和生物学数据，对体长世代分析 (Length-based cohort analysis, LCA) 模型和基于生物量的体长世代分析 (Biomass-based length-cohort analysis, B-LCA) 模型进行性能检验和敏感性分析，并利用蒙特卡洛 (Monte Carlo) 方法估算模型参数、秋刀鱼资源量、捕捞死亡系数以及最大持续产量。结果表明：1) 在5、10和15 mm体长间隔下，LCA和B-LCA模型均表现出优秀的拟合能力，且在5 mm体长间隔下，2种模型的拟合能力均更强；2) LCA模型对于以尾数为单位的渔业数据表现更佳，B-LCA模型对于以质量为单位的渔业数据表现更佳；3) LCA和B-LCA模型对生长因子(b)、渐近体长(L∞)的变化均较敏感，且对b的敏感程度更高；4) LCA模型估算的2014—2018年秋刀鱼平均资源质量约为65.93×104~171.51×104 t，捕捞死亡系数为0.529 2，最大持续产量为37.73×104 t，而B-LCA模型估算的平均资源质量约为47.88×104~126.25×104 t，捕捞死亡系数为0.540 5，最大持续产量为33.02×104 t。2种模型估算的最大持续产量均低于北太平洋渔业委员会 (North Pacific Ocean Commission, NPFC)各成员国年均产量 (40.98×104 t)，表明2014—2018年秋刀鱼资源处于过度捕捞状态。
Abstract:Being one of the primary fishing targets in Chinese pelagic fishing, Pacific saury (Cololabis saira) is distributed in the subtropical to temperate waters of the northwest Pacific Ocean. In order to explore its resource status, according to the catch at size and biological data of Northwest Pacific saury from 2014 to 2018, we conducted a performance test and a sensitivity analysis on length-based cohort analysis (LCA) model and biomass-based length-cohort analysis (B-LCA) model. Besides, we applied Monte Carlo method to estimate the model parameters, resource quantity, fishing mortality coefficients and maximum sustainable yield of Pacific saury. The results show that: 1) The LCA model and B-LCA model exhibited excellent fitting abilities at 5, 10 and 15 mm length intervals, with stronger fitting abilities at 5 mm interval. 2) LCA model performed better for fishery data in units of number, while B-LCA model performed better for fishery data in units of mass. 3) Both LCA and B-LCA models were sensitive to changes in growth factor ($ b $) and asymptote length ($ {L}_{\infty } $), with higher sensitivity to b. 4) The average resource mass of Pacific saury from 2014 to 2018 estimated by LCA model was about 65.93×104−171.51×104 t; the fishing mortality coefficient was 0.529 2; the maximum sustainable yield was 37.73×104 t. The average resource mass estimated by B-LCA model was about 47.88×104−126.25×104 t; the fishing mortality coefficient was 0.540 5; the maximum sustainable yield was $ 33.02\times {10}^{4} $ t. The maximum sustained production estimated by both models was lower than the average annual production of NPFC (North Pacific Ocean Commission) member countries ($ 40.98\times {10}^{4} $ t), indicating that the Pacific saury resources had been overfished from 2014 to 2018.
-
Keywords:
- Cololabis saira /
- Catch at size /
- LCA model /
- B-LCA model /
- Resource quantity /
- Maximum sustainable yield
-
中西太平洋鲣 (Katsuwonus pelamis) 是世界上重要的渔业资源，为全球人口供应动物性蛋白质，并在出口创汇、渔业转型升级及远洋渔业相关产业发展中作出了重要贡献。金枪鱼围网渔业在我国中西太平洋远洋渔业的地位举足轻重，而鲣是围网渔业中的重要捕捞对象[1]。2020年中西太平洋海域捕捞量的72%为金枪鱼围网捕捞，其中鲣占围网总捕捞量的82%[2]。多数研究认为，在影响鲣资源的重要因素中海洋环境因素不可忽视[3-5]，因此，开展有关海洋环境与鲣资源关系的研究具有重要意义。
目前一些研究方法被广泛应用于探讨鲣资源空间分布与海洋环境之间的关系，如随机森林法则[6](Random Forest, RF)、BP神经网络[7-8]，广义加性模型[9] (Generalized Additive Model, GAM) 和最大熵模型[10-11] (Maximum Entropy Model, MaxEnt) 等。然而，这些模型并未考虑海洋环境因子影响的空间异质性问题。针对传统线性回归模型所忽略的空间异质性问题，地理加权回归模型[12-13] (Geographically Weighted Regression, GWR) 在一定程度上有所改进，但该模型只能反映各变量的平均尺度，未考虑不同变量的空间异质性尺度差异，无法体现不同环境因子的多尺度效应，因此也存在一定的估计偏差。
近年来，在分析影响因素空间异质性的实证研究中，逐渐纳入多尺度地理加权回归模型[14] (Multi-scale Geographically Weighted Regression, MGWR)，随着相关的统计推断不断地补充完善，MGWR模型已较普遍地应用于研究中[15-16]，且已证明其在空间尺度差异和异质性研究上有较好的拟合效果[17]。MGWR模型考虑了各影响因子的尺度差异，本研究采用该模型方法，运用到中西太平洋鲣渔获率与海洋环境因子关系研究中，同时为探究MGWR模型的精度，选取了GAM和GWR模型为对照，对比分析这3种模型的拟合优度，探讨了不同环境因子对中西太平洋鲣渔获率空间分布的影响，选出最合适的模型，为鲣资源的养护管理和合理开发利用以及我国金枪鱼围网渔船的生产提供参考依据。
1. 材料与方法
1.1 数据来源
1.1.1 渔业数据
渔业数据来源于中西太平洋渔业委员会 (Western and Central Pacific Fisheries Commission, WCPFC) 所公布的中西太平洋海域围网渔业捕捞数据，本文选取2005—2019年140°E—160°W和15°S—15°N海域范围内的捕捞对象为鲣的数据进行研究，时间分辨率为月，空间分辨率为1°×1°，包括年份、月份、经纬度、捕捞天数及渔获量等。
1.1.2 环境数据
选取与鲣活动和栖息相关的环境因子，包括不同水层 (5、55 、100 、150、200 m) 的温度、盐度、东西向和南北向海水流速，以及净初级生产力。所有环境因子均来自参与第六次国际耦合模式比较计划 (Coupled Model Intercomparison Project, CMIP6) (https://esgf-node.llnl.gov/search/cmip6/) 气候-地球系统模式研发的美国国家大气科学研究中心 (National Center for Atmospheric Research, NCAR) 及美国国家大气海洋局地球流体动力学实验室 (National Oceanic and Atmospheric Administration-Geophysical Fluid Dynamics Laboratory, NOAA-GFDL)，其中2014年及以前的环境数据用于评估CMIP6地球系统模式在历史时期的模拟表现，2015年以后使用了最新的综合评估模型和排放数据[18]， 本文选取其中的2005—2019年的环境数据，并通过Matlab R2018b软件将其与渔业数据进行匹配。时间尺度为月，空间尺度为1°×1°。
1.2 数据预处理
1) 计算2005—2019年1°×1°各单元渔区内的累计渔获量和累计作业时间 (d)，从而获得2005—2019年各单元渔区内渔获率，即名义单位捕捞努力量渔获量 (Catch per unit effort, CPUE)。
2) 考虑到捕捞时的效率容易受各种因素的影响而产生差异，本文采用名义CPUE表示渔获率的分布，其计算公式为：
$$ {Y}_{\mathrm{C}\mathrm{P}\mathrm{U}\mathrm{E}(i,\;j)}=\frac{{{U}_{\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{c}\mathrm{h}}}_{(i,\;j)}}{{{f}_{\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}}}_{(i,\;j)}} $$ (1) 式中：i, j表示作业渔船的经纬度位置；YCPUE表示名义单位捕捞努力量渔获量(t·d−1)；
$ {U}_{\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{c}\mathrm{h}} $ 表示对应地点作业渔船的累计渔获量 (t)；$ {\,f}_{\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}} $ 表示在该作业位置累计作业时间 (d)。1.3 模型与方法
1.3.1 环境因子的选取
当各因子间存在多重共线性时，会出现过拟合而导致模型泛化能力降低。本文采用方差膨胀因子 (Variance inflation factor, VIF) 判断各海洋环境变量是否存在多重共线性，提取VIF<7.5的环境变量[19]。
1.3.2 GAM
GAM是一种全局回归模型，它使用未指定的(即非参数的)平滑函数来代替线性协变量函数。由于该方法可以有效地分析环境因子和渔获率之间的非线性关系[20-21]，因此被广泛运用于渔业研究中。利用GAM模型建立CPUE与环境因子之间的关系式为：
$$ \begin{array}{c} \ln \left(Y_{\mathrm{CPUE}}\right) \sim S(\mathrm{SST})+S(\mathrm{S} 100)+S(\mathrm{U} 55)+S(\mathrm{V} 55)+\\ S(\mathrm{NPP})+ \epsilon _{\mathrm{GAM}} \\[-12pt] \end{array}$$ (2) 式中：S(.)为各环境因子的样条平滑函数；
$ \epsilon $ GAM为模型拟合残差。1.3.3 GWR和MGWR
GWR模型是在传统的全局回归模型的基础上，将数据的地理位置纳入回归参数，在估计局部参数时考虑相邻点的空间权重[22]。其模型表达式如下：
$$ {y}_{i} = {\sum} _{j=1}^{k}{\beta }_{j}\left({u}_{i},{v}_{i}\right){x}_{ij} + {\epsilon }_{i} $$ (3) 与经典GWR模型相比，MGWR模型允许每个变量各自不同的空间平滑水平，这降低了估计的偏误，同时也产生了更真实有用的空间过程模型，MGWR模型的计算公式如下：
$$ y_i={\sum}_{j=1}^k \beta_{\mathrm{bw} j}\left(u_i, v_i\right) x_{i j}+ \epsilon _i$$ (4) 式(3)、式(4)中：
$ {x}_{ij} $ 是变量${x}_{j} $ 在$i $ 点的值；$({u}_{i},{v}_{i}) $ 代表样本点的空间地理位置$ ；{\beta }_{j}\left({u}_{i},{v}_{i}\right) $ 是i点上的第j个回归参数，是地理位置的函数，当j=0时，$ {\beta }_{0}({u}_{i},{v}_{i}) $ 为i点的回归常数；k为回归系数的总个数；$ {\beta }_{\mathrm{b}\mathrm{w}j} $ 代表了第 j 个影响因子回归系数使用的最佳带宽；$ {\epsilon }_{i}$ 为随机误差。MGWR模型的各回归系数$ {\beta }_{\mathrm{b}\mathrm{w}j} $ 均基于局部回归所得，且带宽具备特异性[23]，在经典GWR中，$ {\beta }_{j} $ 所有影响因子的带宽都相同。MGWR模型的核函数和带宽选择准则也与经典GWR模型一致，本文采用Gauss核函数，带宽选取准则为修正的赤池信息准则 (Corrected Akaike InformationCriterion, AICc)。对比可知，MGWR通过推导出响应变量和不同预测变量之间条件关系的单独带宽，允许不同的过程在不同的空间尺度上运行。MGWR使用反向拟合算法进行校准，用GWR参数估计初始化反拟合过程。基于这些初始值，校准过程以迭代的方式工作，在每次迭代中，所有的局部参数估计和最优带宽都被评估。当连续迭代的参数估计的差值收敛于指定阈值时，迭代终止。本研究中，收敛阈值取为10−5。
本文MGWR模型的计算基于美国亚利桑那州立大学空间分析研究中心 (SPARC) 开发的MGWR 2.2软件 (https://sgsup.asu.edu/SPARC)，该软件在进行模型拟合时可选择对环境变量进行标准化处理，其表达式为：
$$ {Z}_{ij}=({X}_{ij}-\overline{X})/{\sigma } $$ (5) 式中：
$ {Z}_{ij} $ 为标准化后的变量值；$ {X}_{ij} $ 为实际变量值；$ \overline{X} $ 和$ \sigma $ 分别为实际数据的均值和标准差。采用ArcGIS 10.5软件制作地图。1.3.4 模型性能评估
为评估模型的性能，本文对比了GAM、GWR和MGWR模型的AICc、残差平方和 (Residual sum of squares, RSS)、拟合优度 (R²) 和校正后的拟合优度 (Adjusted R²)。较低的AICc值、RSS值和较高Adjusted R2值意味着模型具有较好的拟合效果。
2. 结果
2.1 环境因子选取
通过多重共线性诊断，筛选出以下5个海洋环境因子，结果见表1。
表 1 解释变量间方差膨胀因子Table 1. Variance inflation factor among explanatory variables解释变量
Explanatory variable单位
Unit方差膨胀因子
VIF海表面温度
Sea surface temperature, SSTK 1.558 100 m 水深盐度
Sea water salinity at 100 m depth, S100‰ 1.090 净初级生产力
Net primary productivity, NPPmol·(m2·s)−1 1.858 55 m 东西向海水流速
Mean zonal (E-W) current velocity at
55 m depth, U55m·s−1 1.007 55 m 南北向海水流速
Mean meridional (N-S) current
velocity at 55 m depth, V55m·s−1 1.796 2.2 模型性能评估
GAM、GWR和MGWR模型的计算结果参数见表2。与GAM模型相比，GWR模型AICc值略高，调整后的R²显著提升；两个模型的残差平方和 (Residual Sum of Squares, RSS) 对比显示，GWR模型的RSS显著降低，表明GWR模型的拟合度高于GAM模型。在比较GWR模型和MGWR模型的参数时，发现MGWR的RSS和AICc值明显较低，显示出MGWR模型的结果优于GWR模型。此外，MGWR模型调整后的R²从0.846增至0.871，进一步证明了MGWR模型的拟合优度最好。总体而言，MGWR模型的性能明显更好，可以用来解释CPUE的空间分布。
表 2 GAM、GWR 和 MGWR 不同回归模型性能评价对比Table 2. Comparison of statistical parameters of different linear regression models (GAM, GWR and MGWR)参数
Parameter广义加性
模型
GAM地理加权
回归
GWR多尺度地理
加权回归
MGWR残差平方和 RSS 254 400.400 148.607 129.136 修正的赤池信息准则 AICc 1 477.598 1 572.437 1 287.873 拟合优度 R² 0.369 0.879 0.895 校正后的拟合优度
Adjusted R²0.273 0.846 0.871 2.3 模型尺度结果
MGWR模型的结果能直接表现出不同海洋环境因子的差异化作用尺度(表3)，而GWR模型仅能反映各变量作用尺度的平均状态。经典GWR模型的带宽为65，为样本总数量的5.29%。MGWR模型显示各海洋环境因子的作用尺度有很大差异，最佳拟合带宽由小到大依次为NPP、U55、S100、SST、V55。截距表示在其他环境因子确定的情况下，渔获率所受到的影响，其作用尺度即带宽为48，占样本总数量的3.91%，较低于其他变量的作用尺度。SST的作用尺度也较小，为54，占样本总数的4.40%，说明海表面温度存在的空间异质性较大。S100和U55作用尺度均为48，占样本总数3.91%，表明这两个环境因子也具有较大的空间异质性。NPP作用尺度最小，占样本总数量的3.51%，表明空间异质性极强。V55带宽为1 227，属于全局尺度，其系数在空间上表现为平稳缓和，几乎不存在空间异质性。
表 3 MGWR 与 GWR 带宽对比结果Table 3. Bandwidth comparison between classical MGWR and GWR models变量
Variable多尺度地理
加权回归
MGWR占比
Proportion/
%地理加权
回归
GWR占比
Proportion/
%截距 Intercept 48 3.91 65 5.29 海表面温度 SST 54 4.40 65 5.29 100 m 水深盐度 S100 48 3.91 65 5.29 55 m 东西向海水
流速 U5548 3.91 65 5.29 55 m 南北向海水
流速 V551 227 99.84 65 5.29 净初级生产力 NPP 44 3.51 65 5.29 2.4 影响因素的空间异质性
在CPUE与环境影响因素关系的实际讨论中，GWR和MGWR模型都考虑了空间非平稳性特点及空间尺度问题，但MGWR模型为每个影响因素纳入了不同的空间尺度，充分利用每个变量的不同带宽(表3)，对中西太平洋鲣的渔获率进行了更精确的回归分析。MGWR模型的结果中，每个海洋环境因子有其特定的回归系数。各环境因子的具体回归系数信息见表4。
表 4 基于MGWR的各环境因子的局部系数统计描述Table 4. Statistical description of MGWR local coefficient变量
Variable均值
Mean标准差
SD最小值
Min.最大值
Max.正值比例
Positive ratio/%负值比例
Negative ratio/%显著性检验
P截距 Intercept −0.036 0.350 −0.552 0.939 38.975 61.025 1.00 海表面温度 SST 0.162 0.331 −0.555 1.551 66.802 33.198 <0.05 100 m 水深盐度 S100 0.574 0.952 −0.359 4.649 73.556 26.444 <0.05 55 m 东西向海水流速 U55 −0.180 0.504 −2.215 0.997 32.303 67.697 <0.05 55 m 南北向海水流速 V55 0.069 0.109 −0.142 0.553 80.797 19.203 <0.05 净初级生产力 NPP 0.172 0.351 −0.787 1.806 64.849 35.151 <0.05 本文中MGWR模型各环境因子局部回归系数和显著性的空间分布(图1)显示，SST、S100、U55、V55和NPP 5种环境因子存在显著的空间差异。MGWR各因子局部回归系数统计描述见表4。可以发现，SST在大约2/3区域对鲣渔获率有正向影响，1/3部分为负向影响，170°E为界东西两侧分别为负向影响和正向影响(图1-b)，其中在赤道150°E附近回归系数值最大。影响的平均值为0.162，且标准差较小，说明海表面温度每增加1 ℃，渔获率平均增加0.162个单位。由显著性分布可知，170°E以西海域，SST对鲣渔获率正向影响的显著性较强，180°经线处负向影响最为显著。
![]()
图 1 MGWR 模型各环境因子局部回归系数和显著性的空间分布Figure 1. Spatial distribution of local regression coefficients and significance of each environmental factor in MGWR modelS100正向影响区域近73%，平均值为0.574，远大于其他因子，同时标准差也较大，因此变异系数相对较小，表明影响的空间差异较小。从系数空间分布图(图1-c) 可见，总体上S100对鲣渔获率影响的显著性由南至北逐渐增强，其中有两个较为明显的区域，其一为赤道175°E附近，另一显著区域位于赤道西太平洋“暖池”区域；负向影响主要分布在10°S附近。
U55局部回归系数取值介于−2.215~0.997，且变异系数较大，表明U55对渔获率影响的异质性较强，负值比率约68%，因此主要为负向影响。总体上U55对鲣渔获率的影响呈现环形分布状态，且负向影响范围更加广泛。标准差较大为0.504，说明局部回归系数差异也较大。由图1-d可知，赤道170°E附近为负向影响中心，影响程度由此向四周递减；在140°E—155°E、0°—10°S内该系数变为正值，且根据显著性P值分布可知，此处U55对鲣渔获率的正向影响最为明显，但影响范围较小。
V55局部回归系数介于−0.141~0.553，平均值为0.069，标准差为0.109，正值比率超过80%，可见V55在全局上对渔获率产生正向影响，由图1-e可知，V55的负向影响区域较为分散；但从系数绝对值来看，其影响程度相对最小。
NPP主要在65%研究区域内对渔获率产生正向影响，局部回归系数介于−0.787~1.806，平均值为0.172，标准差为0.351。如图1-f，170°E两侧差异明显，在赤道150°E—160°E附近NPP显著性最显著，表明该处与鲣渔获率有较大相关性。170°E以东区域，局部回归系数主要为负值，对渔获率产生负向影响，且显著性较弱。
总体上，各环境因子对鲣渔获率影响的空间异质性程度可由变异系数表示，本研究中空间异质性程度依次为：U55>SST>NPP>S100>V55。由局部回归系数的正负值比例可知，V55和S100正负值相差比例及正值比例均较大，U55的负值比例最大。各环境因子对鲣渔获率影响的空间异质性呈现东西向差异，SST和NPP对渔获率的影响主要受经度位置影响，总体上170°E两侧会呈现不同的分布规律，S100对渔获率的影响南北纬差异较为突出(图1)。同时，局部回归系数的绝对值大小表示各环境因子对鲣渔获率影响程度大小。本研究显示，各环境因子对渔获率的影响程度依次为：S100>U55>NPP>SST>V55。
2.5 模型预测结果比较分析
由GAM、GWR和MGWR模型预测的中西太平洋鲣渔获率标准化后的空间分布图(图2)可见，不同模型的预测范围有所差异，其中GAM模型预测的渔获率高值区介于140°E—178°W、5°S—2°N之间，预测的渔获率低值区域主要分布在研究海域外围 (图2-a)。MGWR模型预测的渔获率高值区见图2-b，分布在A1和A2海域范围内，渔获率低值分布范围主要在B1、B2、B3及B4附近(图2-b)。与GAM模型相比，MGWR模型预测的渔获率高值区范围较小，而低值区范围相对较广。GWR模型预测的渔获率高值分布区域与MGWR模型的结果几乎一致，渔获率低值区主要为C1、C2和C3 (图2-c)。与GWR模型相比，MGWR模型预测的渔获率高值区域范围更集中。总体而言，由3个模型预测的渔获率空间分布与2005—2019年CPUE空间分布(图2-d)可知，与GAM模型相比，MGWR和GWR模型能更好地反映渔获率的分布状况。此外，由于各海洋环境因子具有空间异质性特点，GAM模型无法明确体现其对结果产生的影响，而MGWR和GWR模型均能较好地体现出这一特点。
![]()
图 2 不同模型渔获潜能标准化预测及 2005—2019 年平均 CPUE 空间分布Figure 2. Standardized catch potential prediction and average CPUE spatial distribution of different models from 2005 to 20193. 讨论
3.1 环境因子对鲣渔获率影响的空间异质性
鲣是大洋性中上层高度洄游性鱼类[24]，其资源丰度、集群、资源的分布及洄游与海洋环境息息相关[25]。本研究中MGWR模型结果显示(图1)，SST、S100、U55、V55、NPP 5种环境因子对鲣渔获率的影响均存在显著的空间差异。S100对鲣渔获率主要为正向影响，S100影响下出现两个较为明显的正向区域，主要原因为赤道170°E附近受热带东太平洋信风影响，生成了巨大的涌升流，因此形成了低温、高盐的冷舌区域[26]；S100影响下另一个高值区域和赤道西太平洋“暖池”区域重合度高，“暖池”区域有高温、低盐的特征，当盐度增加时，净初级生产力提升，从而为鲣提供食物资源[27]，导致产生更显著的空间聚集效果(图1-c)。
U55对研究区域内鲣渔获率主要为负向影响(图1-d)。根据洋流分布可知，本研究区域有多条东西向暖流经过，大部分区域都受到了影响，处于洋流流经区域，因此在区域内东西向海流速度对渔获率分布差异产生了较大影响。本研究中U55正向影响明显的区域处于赤道附近，受赤道暖流影响，洋流方向引导鱼群聚集于此[28]。
NPP对渔获率主要为正向影响。浮游动植物的生长通常会受到NPP的影响[29]，NPP值越大的区域，其环境条件对鲣的生存越有利，不同季度冷暖水团的锋面位置会发生变动，NPP通常会随冷暖水团的锋面移动[30]，因此本研究中NPP主要对鲣渔获率存在正向影响。
SST在170°E以西对渔获率有显著的正向影响。已有研究显示，中西太平洋热带海域是一个“暖池-冷舌海洋生态系统”海域[21,31]，在图1-b中，正向影响最为显著的区域基本在中西太平洋的“暖池”区域内，空间分布范围较为广泛，整体资源热度均较好，其整体环境特征是高温和相对较低的叶绿素浓度，并且“暖池”海域靠近陆地及岛屿，因此近岸上升流也提供了较充足的初级生产力[32]。同时在“冷舌”附近渔获率也受到了一定程度的正向影响，主要是受中东太平洋赤道上升流影响，锋面区域也产生了较丰富的初级生产力。金枪鱼鱼群往往在锋面地带获得丰富的食物资源[29,31]，这和本研究的结果基本吻合。
总体上，V55对鲣渔获率影响的空间异质性程度相对最小(表4)，空间异质性的影响范围较小，在正负影响的差异变化上较为平稳缓和(图1-e)；但V55对渔获率的正向影响相对最为明显，主要体现在靠近巴布亚新几内亚的海域及180°附近(图1-e)。U55空间异质性程度最大，但目前较少有研究使用海水流速等动力学因子进行分析，未来的研究建议考虑相关的海洋动力学因子。SST和NPP也有较明显的影响。S100对渔获率的影响在纬度变化上也有东西向差异，但南北纬向差异更为突出，NPP和SST对渔获率的影响主要受经度位置影响，总体上170°E两侧会呈现不同的分布规律。这是由于170°E以东海域受热带东太平洋信风影响，从而产生了涌升流，形成了温度低、盐度高、初级生产力高的冷舌区。暖池边缘生成的由底层海水上泛且带有大量营养物质的辐合区，聚集了大量的微型浮游动物及浮游植物，使鲣资源拥有了丰富的饵料[33]，形成了良好的产卵场和索饵场。在170°E以西鲣的空间聚集区附近海域，海表面温度相对较高，西部暖锋有两个“暖池”向东凸出，有利于锋面的形成，产生的丰富食物资源为鲣提供了良好的生存条件，鲣的资源量因此较为丰富。
3.2 环境因子对鲣渔获率影响程度
本研究中MGWR模型各环境因子回归系数的绝对值显示，各环境因子对渔获率影响程度大小不同。5个影响因子中，S100对渔获率的影响程度最大，其次为U55、NPP和SST，V55的影响程度相对较小。盐度对鲣渔场影响的研究表明，盐度在海洋变化中充当重要的物理参数[34-35]。大洋环流中海水盐度是有效的示踪物，且通过对海水密度的调节变化，在海洋热量的垂直输送上产生了重要影响[36]，海表面温度异常也间接受到影响[37]。本研究中热带中西太平洋海域水分和温度条件充足，盐分对表层营养物质至关重要，从而对鲣资源产生重要影响。同时混合层深度受温度和盐度作用，对鲣资源也产生了一定影响[38]。U55对渔获率大小的影响程度也相对较高。Meehl[28]认为，东西向海水流速在热带地区有最大的季节性变化，海流的变化会引起盐度、溶解氧、温度等因素发生改变，使诱饵产生聚集，海流方向可以引导鱼类游向诱饵来源，鱼群更容易定位诱饵。受信风的影响，东太平洋产生了强烈的涌升流，形成了盐度高、温度低等特征的冷舌区域。冷暖流相遇形成的交汇海域浮游动植物较为充足，是鲣理想的索饵场。因此，本研究认为NPP也是影响鲣渔获率的重要因子之一。张小龙等[39]认为，由于148°E—165°E、6°N—6°S附近海域存在两个向东凸出的“暖池”易形成温度锋面，为鲣的生存提供了丰富的饵料资源，形成了明显的空间聚集现象，因此SST在一定程度上也影响了中西太平洋鲣渔获率。以往研究认为SST是影响鲣渔获率的重要因子，但本研究采用平均状态下的环境和资源数据，因此SST的影响程度较次于S100和NPP。同时从空间形态观察来看，S100、NPP和SST这三者区域差异的影响模式相似。
3.3 环境因子多尺度效应
本研究中，由于MGWR模型考虑了环境因子的多尺度效应，结果发现NPP作用尺度较小，其次为S100和U55，SST的作用尺度较大，V55最大，为全局尺度。NPP是主要表现食物的因子，其作用尺度较小可能是由于鲣资源捕食范围受该区域冷暖水团的影响，而冷暖水团边界范围较狭窄，所以其影响的空间范围也相对较小。鲣是喜温性鱼类，主要分布在水温高于29.5 ℃的区域内[34]，本研究区域内暖水团范围较广，因此SST的作用尺度较大。近岸上升流和盐度锋面均对营养盐有一定影响[32]，但由于位置较为固定，因此S100的作用尺度也受到了一定的空间范围限制。
总体上，SST、NPP对鲣渔获率空间分布的影响相似，受到经度位置的影响，170°E两侧呈现不同的规律，其中170°E以西主要为正向影响，且具有明显的岛屿属性；170°E以东主要为负向影响。这是由于170°E以东海域在热带东太平洋信风的作用下，形成温度低、盐度高的冷舌区域；170°E以西海域，正向影响最为显著区域基本在中西太平洋的“暖池”区域内，整体资源热度均较好，整体上温度较高、叶绿素浓度相对较低，并且该区域靠近陆地和岛屿，因此近岸上升流也为其提供了丰富的营养物质。S100南北向差异较为明显，而V55没有明显的东西或南北向差异。U55正向影响区域较小，主要集中在热带岛屿附近，原因可能是海流的变化引起岛屿附近盐度、溶解氧、温度等因素发生变化，使饵料产生了聚集[28]。
4. 小结与展望
目前，有关环境与鲣关系的研究主要采用MaxEnt模型[12]、GAM[40]模型等，但这些模型均未考虑空间异质性问题。GWR模型可以通过结合不同地理位置的空间属性差别，来探索不同环境因子在时空上的差异性影响，并获得海洋环境与资源间的时空异质性特征。本研究在GWR模型基础上引入MGWR模型，加入了多尺度概念，能捕捉到不同变量的不同影响尺度，考虑了不同环境因子在影响尺度上的差异性，从而避免了捕获太多的噪声和偏误。因此，是否考虑影响因素的空间尺度会对模型的结果和分析产生巨大影响。由于本研究采用平均状态的环境和资源数据，所用的MGWR模型未顾及时间序列对海洋资源的影响，建议未来的研究在此基础上加入时空地理加权回归模型，综合时间和空间两个维度进行分析。
-
![]()
图 3 基于尾数模拟数据LCA和B-LCA模型的估算值与模拟值对比图
Figure 3. Comparison between estimated and simulated values of LCA and B-LCA models based on number-based simulated data
![]()
图 4 基于质量模拟数据LCA和B-LCA模型的估算值与模拟值对比图
Figure 4. Comparison between estimated and simulated values of LCA and B-LCA models based on mass-based simulated data
![]()
图 5 基于LCA和B-LCA模型估算的平均资源质量与捕捞死亡系数在体长组间的变化趋势
Figure 5. Trends of average resource mass and fishing mortality coefficient between body length class estimated based on LCA and B-LCA models
表 1 秋刀鱼生物学参数
Table 1 Biological parameters of Pacific saury
参数
Parameter符号
Symbol数值
Value条件因子 Condition factor a 2.468 4×10−6 生长因子 Growth factor b 3.110 7 自然死亡系数 Coefficient of natural mortality M 0.328 9 捕捞死亡系数 Coefficient of fishing mortality F 1.271 1 总死亡系数 Total mortality coefficient Z 1.600 0 渐进体长 Asymptote length L∞ 360.230 0 生长参数 Growth parameter K 0.360 0 
下载: 导出CSV
表 2 模型性能指标
Table 2 Model performance indicators
间隔
Interval/mm模拟数据
Simulated data模型
Model性能指标
Model efficiency5 基于资源尾数
Number-basedLCA 1.000 0 B-LCA 0.995 9 基于资源质量
Mass-basedLCA 0.995 5 B-LCA 1.000 0 10 基于资源尾数
Number-basedLCA 1.000 0 B-LCA 0.986 4 基于资源质量
Mass-basedLCA 0.983 3 B-LCA 1.000 0 15 基于资源尾数
Number-basedLCA 0.999 9 B-LCA 0.972 6 基于资源质量
Mass-basedLCA 0.962 7 B-LCA 1.000 0
下载: 导出CSV
表 3 敏感性分析结果
Table 3 Sensitivity analysis results
模拟数据
Simulated data参数
Parameter基准方案
Base case变化
Change敏感性方案
Sensitivity case敏感性指数
Sensitivity index$ {W}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}^{\mathrm{L}\mathrm{C}\mathrm{A}} $ $ {W}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}^{\mathrm{B}{\text{-}}\mathrm{L}\mathrm{C}\mathrm{A}} $ $ {W}_{\mathrm{s}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y}}^{\mathrm{L}\mathrm{C}\mathrm{A}} $ ${W}_{\mathrm{s}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{i}\mathrm{t}\mathrm{y} }^{\mathrm{B}{\text{-}}\mathrm{L}\mathrm{C}\mathrm{A} }$ LCA B-LCA 基于资源尾数 a 1.604 7 1.487 2 +10% 1.765 2 1.635 9 1.000 0 0.194 6 Number-based 1.604 7 1.487 2 −10% 1.444 2 1.338 5 1.000 0 1.658 9 基于资源质量 0.088 4 0.081 9 +10% 0.088 4 0.081 9 0.000 0 −0.732 2 Mass-based 0.088 4 0.081 9 −10% 0.088 4 0.081 9 0.000 0 0.732 2 基于资源尾数 b 1.604 7 1.487 2 +10% 8.448 0 7.775 2 42.645 3 38.452 7 Number-based 1.604 7 1.487 2 −10% 0.305 6 0.285 2 8.095 7 8.222 6 基于资源质量 0.088 4 0.081 9 +10% 0.095 7 0.088 1 0.831 5 −0.031 2 Mass-based 0.088 4 0.081 9 −10% 0.081 8 0.076 4 0.744 3 1.361 1 基于资源尾数 $ K $ 1.604 7 1.487 2 +10% 1.786 5 1.657 7 1.133 1 0.330 0 Number-based 1.604 7 1.487 2 −10% 1.786 5 1.657 7 1.116 8 1.777 6 基于资源质量 0.088 4 0.081 9 +10% 0.098 4 0.091 3 1.135 2 0.331 8 Mass-based 0.088 4 0.081 9 −10% 0.098 4 0.091 3 −1.135 2 −0.331 8 基于资源尾数 $ {L}_{{\infty } } $ 1.604 7 1.487 2 +10% 2.080 8 1.935 1 2.966 7 2.058 7 Number-based 1.604 7 1.487 2 −10% NA NA NA NA 基于资源质量 0.088 4 0.081 9 +10% 0.114 7 0.106 7 2.978 2 2.069 3 Mass-based 0.088 4 0.081 9 −10% NA NA NA NA 基于资源尾数 $ M $ 1.604 7 1.487 2 +10% 1.605 3 1.487 7 0.003 5 −0.728 9 Number-based 1.604 7 1.487 2 −10% 1.604 2 1.486 7 0.003 4 0.735 4 基于资源质量 0.088 4 0.081 9 +10% 0.088 4 0.081 9 0.003 3 −0.729 2 Mass-based 0.088 4 0.081 9 −10% 0.088 4 0.081 9 0.003 1 0.735 1 基于资源尾数 $ F $ 1.604 7 1.487 2 +10% 1.764 7 1.635 5 0.997 0 0.191 8 Number-based 1.604 7 1.487 2 −10% 1.444 7 1.338 9 0.997 0 1.656 1 基于资源质量 0.088 4 0.081 9 +10% 0.081 1 0.075 1 −0.827 6 −1.507 9 Mass-based 0.088 4 0.081 9 −10% 0.097 1 0.090 1 −0.981 8 −0.192 1 注：NA表示缺失值。 Note: NA indicates missing values.
下载: 导出CSV
表 4 模型参数的后验分布
Table 4 Posterior distribution of model parameters
参数
Parameter模型
Model95%置信区间
95% confidence interval均值
Mean标准差
SD潜在尺度缩减因子
Psrf下限 Lower 上限 Uppe $ M $ LCA 0.264 13 0.392 13 0.328 93 0.032 697 1.000 1 B-LCA 0.265 26 0.394 98 0.328 75 0.033 046 1.000 0 $ F $ LCA 1.018 6 1.520 7 1.270 8 0.128 04 1.000 1 B-LCA 1.018 1 1.516 5 1.270 5 0.127 10 1.000 2 $ K $ LCA 0.292 02 0.432 27 0.360 04 0.035 822 1.000 0 B-LCA 0.288 91 0.430 53 0.360 27 0.036 064 1.000 2 $ {L}_{\infty } $ LCA 290.67 429.86 360.24 35.853 1.000 1 B-LCA 288.41 429.55 360.52 36.011 1.000 1 $ b $ LCA 2.504 2 3.729 1 3.110 8 0.311 77 1.000 1 B-LCA 2.505 1 3.722 3 3.111 2 0.310 42 1.000 0 $ a $ LCA 1.99×10−6 2.96×10−6 2.47×10−6 2.47×10−7 1.000 2 B-LCA 1.98×10−6 2.94×10−6 2.47×10−6 2.48×10−7 0.999 99
下载: 导出CSV
表 5 秋刀鱼平均资源质量、捕捞死亡系数及最大持续产量
Table 5 Average resource mass, fishing mortality coefficient and maximum sustainable yield of Pacific saury
年份
Year总产量
Wtc/104 tLCA B-LCA F W/104 t W/104 t F W/104 t WMSY/104 t 2014 63.02 0.422 6 171.51 59.72 0.427 5 126.25 52.26 2015 35.89 0.459 8 86.51 32.17 0.466 4 61.81 28.11 2016 36.17 0.458 7 86.35 32.29 0.462 5 61.68 28.22 2017 26.26 0.667 5 65.93 23.97 0.687 8 47.88 21.00 2018 43.59 0.637 6 113.86 40.52 0.658 1 83.29 35.49 均值 Mean 40.99 0.529 2 104.83 37.73 0.540 5 76.18 33.02
下载: 导出CSV
-
[1] IWAHASHI M, ISODA Y, ITO S I, et al. Estimation of seasonal spawning ground locations and ambient sea surface temperatures for eggs and larvae of Pacific saury (Cololabis saira) in the western North Pacific[J]. Fish Oceanogr, 2006, 15(2): 125-138. doi: 10.1111/j.1365-2419.2005.00384.x
[2] 王茜, 崔雪森. 国际渔业动态[J]. 渔业信息与战略, 2020, 35(4): 327-332. [3] 史登福, 张魁, 陈作志. 基于生活史特征的数据有限条件下渔业资源评估方法比较[J]. 中国水产科学, 2020, 27(1): 12-24. [4] COSTELLO C, OVANDO D, HILBORN R, et al. Status and solutions for the world's unassessed fisheries[J]. Science, 2012, 338(6106): 517-520. doi: 10.1126/science.1223389
[5] HORDYK A, ONO K, VALENCIA S, et al. A novel length-based empirical estimation method of spawning potential ratio (SPR), and tests of its performance, for small-scale, data-poor fisheries[J]. ICES J Mar Sci, 2015, 72(1): 217-231. doi: 10.1093/icesjms/fsu004
[6] RUDD M B, THORSON J T. Accounting for variable recruitment and fishing mortality in length-based stock assessments for data-limited fisheries[J]. Can J Fish Aquat Sci, 2018, 75(7): 1019-1035. doi: 10.1139/cjfas-2017-0143
[7] WANG J T, YU W, CHEN X J, et al. Stock assessment for the western winter-spring cohort of neon flying squid (Ommastrephes bartramii) using environmentally dependent surplus production models[J]. Sci Mar, 2016, 80(1): 69-78.
[8] TRUESDELL S B, BENCE J R, SYSLO J M, et al. Estimating multinomial effective sample size in catch-at-age and catch-at-size models[J]. Fish Res, 2017, 192: 66-83. doi: 10.1016/j.fishres.2016.11.003
[9] JONES R. Assessing the long-term effects of changes in fishing effort and mesh size from length composition data[J]. ICES J Mar Sci, 1974, 33: 1-13.
[10] 詹秉义. 渔业资源评估[M]. 北京: 中国农业出版, 1995: 312-314. [11] 周永东, 徐汉祥. 应用体长股分析法估算东海海鳗资源量[J]. 浙江海洋学院学报(自然科学版), 2007, 26(4): 399-403. [12] 周永东, 张洪亮, 徐汉祥, 等. 应用体长股分析法估算东海区日本鲭资源量[J]. 浙江海洋学院学报(自然科学版), 2011, 30(2): 91-94. [13] 吴斌, 方春林, 贺刚, 等. 运用体长股法初步估算湖口江段短颌鲚资源量[J]. 湖北农业科学, 2014, 53(16): 3866-3869. [14] ZHANG C I, SULLIVAN P J. Biomass-based cohort analysis that incorporates growth[J]. Trans Am Fish Soc, 1988, 117(2): 180-189. doi: 10.1577/1548-8659(1988)117<0180:BBCATI>2.3.CO;2
[15] ZHANG C I, MEGREY B A. A simple biomass-based length-cohort analysis for estimating biomass and fishing mortality[J]. Trans Am Fish Soc, 2010, 139(3): 911-924. doi: 10.1577/T09-041.1
[16] 朱清澄, 花传祥. 西北太平洋秋刀鱼渔业[M]. 北京: 海洋出版社, 2017: 1-2. [17] TAMURA T, FUJISE Y. Geographical and seasonal changes of the prey species of minke whale in the Northwestern Pacific[J]. ICES J Mar Sci, 2002, 59(3): 516-528. doi: 10.1006/jmsc.2002.1199
[18] SHIMIZU Y, TAKAHASHI K, ITO S I, et al. Transport of subarctic large copepods from the Oyashio area to the mixed water region by the coastal Oyashio intrusion[J]. Fish Oceanogr, 2009, 18(5): 312-327. doi: 10.1111/j.1365-2419.2009.00513.x
[19] SMITH A D, BROWN C J, BULMAN C M, et al. Impacts of fishing low-trophic level species on marine ecosystems[J]. Science, 2011, 333(6046): 1147-1150. doi: 10.1126/science.1209395
[20] PIKITCH E K, ROUNTOS K J, ESSINGTON T E, et al. The global contribution of forage fish to marine fisheries and ecosystems[J]. Fish Fish, 2014, 15(1): 43-64. doi: 10.1111/faf.12004
[21] North Pacific Fisheries Commission. 7th Meeting Report[R]. Tokyo: NPFC, 2022.
[22] NAKAYA M, MORIOKA T, FUKUNAGA K, et al. Growth and maturation of Pacific saury Cololabis saira under laboratory conditions[J]. Fish Sci, 2010, 76: 45-53. doi: 10.1007/s12562-009-0179-9
[23] BAITALIUK A, ORLOV A, ERMAKOV Y K. Characteristic features of ecology of the Pacific saury Cololabis saira (Scomberesocidae, Beloniformes) in open waters and in the northeast Pacific Ocean[J]. J Ichthyol, 2013, 53(11): 899-913. doi: 10.1134/S0032945213110027
[24] 朱清澄, 杨明树, 高玉珍, 等. 西北太平洋秋刀鱼耳石生长与性成熟度、个体大小的关系[J]. 上海海洋大学学报, 2017, 26(2): 263-270. [25] KIMURA N, OKADA Y, MAHAPATRA K. Relationship between saury fishing ground and sea surface oceanographic features determined from satellite data along the northeastern coast of Japan[J]. J Mar Sci Technol, 2004, 2(2): 1-12.
[26] 张孝民, 石永闯, 李凡, 等. 基于MAXENT模型预测西北太平洋秋刀鱼潜在渔场[J]. 上海海洋大学学报, 2020, 29(2): 280-286. [27] 张培超, 张孝民, 陈丙见, 等. 北太平洋秋刀鱼渔场与水温垂直结构的关系[J]. 河北渔业, 2022(6): 35-39, 44. [28] 谢斌, 汪金涛, 陈新军, 等. 西北太平洋秋刀鱼资源丰度预报模型构建比较[J]. 广东海洋大学学报, 2015, 35(6): 58-63. [29] 朱文涛, 陈新军, 汪金涛, 等. 基于灰色系统的西北太平洋秋刀鱼资源丰度预测[J]. 广东海洋大学学报, 2018, 38(6): 13-17. [30] 孟令文, 朱清澄, 花传祥, 等. 栖息地指数模型在北太公海秋刀鱼渔情预报中的应用[J]. 海洋湖沼通报, 2018(6): 142-149. doi: 10.13984/j.cnki.cn37-1141.2018.06.019 [31] 许巍, 朱清澄, 张先存, 等. 西北太平洋秋刀鱼舷提网捕捞技术[J]. 齐鲁渔业, 2005, 22(10): 43-45. [32] 石永闯, 朱清澄, 张衍栋, 等. 基于模型试验的秋刀鱼舷提网纲索张力性能研究[J]. 中国水产科学, 2016, 23(3): 704-712. [33] PAULY D. On the interrelationships between natural mortality, growth parameters, and mean environmental temperature in 175 fish stocks[J]. ICES J Mar Sci, 1980, 39(2): 175-192. doi: 10.1093/icesjms/39.2.175
[34] ALLEN M S, WALTERS C J, MYERS R. Temporal trends in largemouth bass mortality, with fishery implications[J]. N Am J Fish Manag, 2008, 28(2): 418-427. doi: 10.1577/M06-264.1
[35] HUGHES S E. Stock composition, growth, mortality, and availability of Pacific saury, Cololabis saira, of the northeastern Pacific Ocean[J]. Fish Bull, 1974, 72(1): 121-131.
[36] GUPTA H V, KLING H, YILMAZ K K, et al. Decomposition of the mean squared error and NSE performance criteria: implications for improving hydrological modelling[J]. J Hydrol, 2009, 377(1/2): 80-91.
[37] MORIASI D N, GITAU M W, PAI N, et al. Hydrologic and water quality models: performance measures and evaluation criteria[J]. Transactions of ASABE, 2015, 58(6): 1763-1785. doi: 10.13031/trans.58.10715
[38] SUN X, NEWHAM L T, CROKE B F, et al. Three complementary methods for sensitivity analysis of a water quality model[J]. Environ Model Softw, 2012, 37: 19-29. doi: 10.1016/j.envsoft.2012.04.010
[39] PLUMMER M, BEST N, COWLES K, et al. CODA: convergence diagnosis and output analysis for MCMC[J]. R News, 2006, 6(1): 7-11.
[40] 官文江, 吴佳文, 曹友华. 利用后向预报方法分析印度洋黄鳍金枪鱼资源评估模型[J]. 中国海洋大学学报 (自然科学版), 2020, 50(2): 52-59. [41] KIM T K. T Test as a parametric statistic[J]. Korean J Anesthesiol, 2015, 68(6): 540-546. doi: 10.4097/kjae.2015.68.6.540
[42] RIVARD D. Effects of systematic, analytical, and sampling errors on catch estimates: a sensitivity analysis[J]. Can Spec Publ J Fish Aquat Sci, 1983, 66: 114-129.
[43] LAI H L, GALLUCCI V F. Effects of parameter variability on length-cohort analysis[J]. ICES J Mar Sci, 1988, 45(1): 82-92. doi: 10.1093/icesjms/45.1.82
[44] 李忠炉, 金显仕, 单秀娟, 等. 小黄鱼体长-体质量关系和肥满度的年际变化[J]. 中国水产科学, 2011, 18(3): 602-610. [45] 吴斌, 方春林, 贺刚, 等. FiSAT II 软件支持下的体长股分析法探讨[J]. 南方水产科学, 2013, 9(4): 94-98. [46] HOLLOWED A B, IANELLI J N, LIVINGSTON P A. Including predation mortality in stock assessments: a case study for Gulf of Alaska walleye pollock[J]. ICES J Mar Sci, 2000, 57(2): 279-293. doi: 10.1006/jmsc.1999.0637
计量
- 文章访问数: 175
- HTML全文浏览量: 98
- PDF下载量: 53
粤公网安备 44010502001741号
