基于贝叶斯状态空间产量模型的大西洋黄鳍金枪鱼资源评估

田志盼，田思泉,2,3，戴黎斌，麻秋云,2,3*

(1.上海海洋大学 海洋科学学院，上海 201306；2.国家远洋渔业工程技术研究中心，上海 201306；3.大洋渔业资源可持续开发教育部重点实验室，上海 201306)

1 引言

黄鳍金枪鱼（Thunnus albacares）是高度洄游类的鱼种，广泛分布于三大洋的热带和亚热带海域，资源量相对丰富且价值较高，是金枪鱼渔业中重要的经济鱼种。随着人类需求的逐渐上升，其受到的捕捞威胁也不断增加，因此进行科学的资源评估进而制定合理的养护管理措施，是实现渔业可持续开发的基础。

目前，国际上对公海金枪鱼实施管辖的是各类区域性渔业管理组织（Regional Fisheries Management Organizations,RFMOs)，大西洋黄鳍金枪鱼由国际大西洋金枪鱼养护委员会（International Commission for Conservation of Atlantic Tunas,ICCAT）进行管理。ICCAT 在对黄鳍金枪鱼的资源评估中，主要使用的评估模型有非平衡剩余产量模型（A Stock Production Model Incorporating Covariates,ASPIC）[1]、年龄结构产量模型（Age-Structured Production Model,ASPM）[2]、实际种群分析（Virtual Population Analysis,VPA）[3]和资源综合模型（Stock Synthesis III,SS3）[4]等，每个模型评估的资源状态不尽相同，ICCAT 综合各模型认为，当前资源处于资源型过度捕捞而无捕捞型过度捕捞状态。其中，剩余产量模型，相较其他模型来说，对渔业数据需求较低（仅需渔获量和种群丰度指数），且能得到最大可持续产量（Maximum Sustainable Yield,MSY）等参考点信息，故是各类渔业资源评估中使用最广泛的模型之一[5]。

除了观测误差，渔业资源种群动态中还存在环境变化等因素产生的过程误差，而ICCAT 当前使用的ASPIC 等产量模型无法估计过程误差，由该类误差产生的不确定性难以被考虑在内。JABBA（Just Another Bayesian Biomass Assessment）是一种基于贝叶斯方法的状态空间产量模型[6]，其中贝叶斯框架可以通过合理的先验信息来降低模型中的不确定性[7-8]，状态空间建模则可以同时估计过程误差和观测误差[9-11]。为此，本文尝试采用JABBA 模型来评估大西洋黄鳍金枪鱼的资源状况，研究状态空间建模对该资源进行评估的适用性，以期为该重要鱼种的科学研究和渔业管理提供更多基础资料和参考信息。

2 材料与方法

2.1 数据来源

本文利用的1950−2017 年渔获量数据来自ICCAT 数据库，渔获量在1990 年达到最高的19.36×104t，2017 年为13.53×104t（图1）。单位捕捞努力量渔获量（Catch Per Unit Effort，CPUE）数据来自ICCAT 黄鳍金枪鱼资源评估会议报告和CPUE 标准化研究报告[12-16]，共计8 个延绳钓船队的标准化CPUE 数据（表1）。虽然各船队根据各自渔业和数据情况采用了不同的CPUE 标准化方法（日本、委内瑞拉和美国为广义线性模型，而乌拉圭1、乌拉圭2、巴西、中国台北1 和中国台北2 则为广义线性混合模型），但应ICCAT 相关工作组的建议，本文资源评估模型在建模时纳入了所有8 个CPUE 数据。

2.2 JABBA 模型

JABBA（版本为v1.1[6]）中运行的Pella-Tomlinson 剩余产量函数形式如下：

式中，SP为剩余产量；r为种群的内禀增长率；K为平衡状态时的未开发资源生物量（即环境容纳量）；B为资源量；m为形状参数。

获得MSY 时的生物量BMSY和捕捞死亡率FMSY分别为

图1 大西洋黄鳍金枪鱼1950−2017 年的年渔获量Fig.1 The annual catch of Atlantic yellowfin tuna from 1950 to 2017

表1 大西洋黄鳍金枪鱼各延绳钓船队标准化CPUE 数据Table 1 Standardized CPUE for each longline fleet of Atlantic yellowfin tuna

捕捞死亡率F定义为

MSY 定义为

B/BMSY＜1 表示当前种群已发生资源型过度捕捞，F/FMSY＞1 表示种群正遭受捕捞型过度捕捞。如果m=2，则SP函数为Scheafer 形式；如果m趋近于1，则为Fox 形式。在JABBA 中，采用默认设置值BMSY/K=0.4，由此得出m=1.2。

过程方程定义如下：

式中，y为年份，Py为y年B与K的比值；ηy为过程误差，且为过程方差，服从逆伽马分布（inverse-gamma (4,0.01)）；Cf,y为y年船队f的渔获量。

JABBA 中观测方程定义如下：

式中，qi为丰度指数i的可捕性系数；εy,i为观测误差，且为观测方差，包含固定项和预测项，预测项服从无信息的逆伽马分布（inverse-gamma(0.001,0.001)）。

本研究中各参数的先验分布设置如下：BMSY/K=0.4；σfix=0.2；初始资源消耗率B1950/K服从对数正态分布，其中值和变异系数分别为1.0 和0.1；r和K的先验信息参考Matsumoto 等[1]的研究结果：假设r服从0.14～0.34 的均匀分布，K服从139.2×104～265.8×104的均匀分布；可捕性系数q为无信息均匀分布。

因Schaefer 的对称形式不符合黄鳍金枪鱼种群动态变化情况[17]，本文只考虑选择Fox 和Pella-Tomlinson 函数。根据不同的CPUE 数据和剩余产量函数，预实验共设置了S1−S8 共8 种方案进行分析（表2）。当均方根误差（Root Mean Squared Error，RMSE）或偏差信息准则（Deviation Information Criteria，DIC）较小时，说明模型拟合效果较好[18]。

当选择所有CPUE 数据并使用Pella-Tomlinson 函数时，得到S1；在S1 预实验基础上，去掉拟合效果差的CPUE 数据得到了S2；在S2 基础上，考虑到ICCAT在CPUE 数据方面的建议[12]—因CPUE 标准化当中未考虑目标鱼种的变化，认为日本延绳钓标准化CPUE数据应该从1976 年开始，舍弃之前的数据得到S3；S3 的CPUE 数据不变，剩余产量函数选择Fox 得到S4；在S4 基础上，考虑对CPUE 数据的敏感性，依次去掉1 条CPUE 数据后得到S5−S8 共4 种方案。

表2 大西洋黄鳍金枪鱼JABBA 模型S1−S8 方案设置Table 2 Different scenarios (S1−S8) of Atlantic yellowfin tuna in JABBA

2.3 回溯性分析

随着渔业数据逐年增加到资源评估中，模型估算结果可能因为出现系统性偏差而导致持续高估或低估的问题称为回溯性问题（Retrospective Problem，RP）。RP 误差的强度主要由Mohn[19]定义的ρ来衡量：

式中，y1、y2分别为整个数据的起始年和结束年，y1：y表示利用y1到y年的数据进行模型估计；X为某一估计的模型参数（如资源生物量或捕捞努力量等）。

如果ρ趋于0，则表明不存在RP；ρ大于0，则存在正RP，即同一年某参数短时间序列的估计值大于整个时间序列的估计值，反之则为负RP[20]。

2.4 敏感性分析

本文通过敏感性分析，研究了种群关键参数K和r的先验分布以及渔获量数据的误报比例对评估结果的影响，进而探讨模型的稳健性。本文分别研究K和r无信息的先验分布和有信息的先验分布（表3）。20 世纪90 年代中后期ICCAT 数据收集上报过程才更加规范[21-22]，这段时间前后数据的可信度存在差异。鉴于此，本文假设1950−1994 年间渔获量数据存在不同程度的误报问题，即上报渔获量占实际渔获量的比例分别设为70%、80%、90%、110%、120%和130%共6 种情况。

2.5 预测分析

ICCAT 在2016 年大西洋黄鳍金枪鱼资源评估会议中，预测分析显示，渔获量低于12×104t 时能使种群到2024 年一直保持健康状态，所以将其总允许可捕量（Total Allowable Catch，TAC）设定为11×104t[22]。因此本研究以11×104t 为基础，设置8.80×104t (80%)、9.35×104t (85%)、9.90×104t (90%)、10.45×104t (95%)、11.00×104t (100%)、12.10×104t (110%)和13.20×104t(120%)共7 个TAC 指标，假设2018 年渔获量为2015−2017 年的平均值，并以2019 年为起始年，预测2019−2027 年的种群动态变化，并以生物量B＞BMSY及种群处于健康状态等的概率评价TAC 指标的管理效果。

3 结果

3.1 模型结果

各方案下得到了JABBA 模型的CPUE 指数趋势（图2）和拟合优度（表4）。在S1 方案下，URU1、URU2、BR 和TAI2 的拟合效果非常差（图2a），存在许多异常值且RMSE 极高（表4），而去除异常值后，拟合效果有较大改善，RMSE 大幅降低（图2b）。S3使用JAP_RE 后，拟合效果略有改善，S4 改用Fox 函数后RMSE 基本不变但DIC 降低。S5、S6、S8 下的拟合效果稍有提升，S7 则变差（表4）。综上所述，且鉴于S4 方案涵盖了更多船队的CPUE 信息，本文将S4 方案作为基础模型来提供资源评估结果。

表3 大西洋黄鳍金枪鱼JABBA 模型中K 和r 有信息和无信息的先验分布设定以及后验分布Table 3 The informative and non-informative prior and posterior distributions for K and r in the JABBA for Atlantic yellowfin tuna

基础模型所有参数后验分布均左右对称且在合理的范围内，说明模型收敛并得到了可靠的结果（图3）。本文求得大西洋黄鳍金枪鱼的MSY 为13.7×104t，BMSY为65.2×104t，B2017略高于BMSY，F2017略低于FMSY（表5）。

3.2 资源状态

随着渔业开发程度的增加，种群由初始（1950 年）的健康状态逐渐进入资源捕捞过度的状态（1997 年前后），随后逐渐恢复，2017 年资源有65%的概率既没有处于资源型过度捕捞状态，也没有处于捕捞型过度捕捞状态，资源状态健康（图4）。20 世纪90 年代和21 世纪初期，种群处于过度捕捞状态（相对捕捞死亡率F/FMSY＞1，而相对生物量B/BMSY＜1（图5））。

以2019 年为起始年，在7 个不同的TAC 目标下，预测分析显示，2019−2027 年资源量均保持增长的趋势（图6）。当TAC 为8.8×104t 时，生物量增长最快，随着TAC 变大，资源量增长速度放缓。风险分析结果显示（表6至表8），TAC为11×104t时，2024年B＞BMSY和资源健康的概率均为84.6%，F＞FMSY的概率为6.2%；当TAC 为13.2×104t 时，2024 年B＞BMSY的概率降低到了69.2%（表6），但F＞FMSY的概率明显增大（28.9%，表7）。

敏感性分析结果表明，当K为无信息先验时，K的估计值略有减小，r估计值略有增大；当r无信息先验时，K的估计值略有增大，r的估计值略有减小（表3）。F2017随渔获量少报程度增大而减小，B2017则与之相反；F2017随多报程度增大而减小，B2017则与之相反（表9）。但不同误报率情况下，B2017/BMSY都大于1，F2017/FMSY都小于1，且资源处于健康状态的概率变化较小。回溯性分析结果表明，当数据逐年减少至2012 年，B、B/BMSY估计值略有减小，F、F/FMSY估计值略有增大，但差别极小（图7）。计算得到B、B/BMSY、F、F/FMSY的ρ值分别为−0.360、0.448、−0.183、0.296。

4 讨论

本文通过贝叶斯状态空间的剩余产量模型，在JABBA 中评估了1950−2017 年大西洋黄鳍金枪鱼的资源状况。当前种群处于没有过度捕捞的健康状态，在ICCAT 当前的TAC 养护管理措施下，2024 年能够达成其保持种群健康状态的养护管理目标。研究结果表明，当使用美国、委内瑞拉、日本（去除1976 年以前）、中国台北1993−2014 年4 个CPUE 数据及Fox函数时，JABBA 模型的拟合效果最佳，评估结果对参数K、r的先验分布和1950−1994 年间的渔获量误报不太敏感，且模型不存在明显的回溯性误差。

1994 年ICCAT 成立工作组对大西洋黄鳍金枪鱼进行评估[22]，之后分别在2000 年、2003 年、2008 年、2011 年和2016 年都进行了资源评估。在2016 年的资源评估中，ASPIC 模型评估认为2014 年大西洋黄鳍金枪鱼处于资源型过度捕捞状态，但没有遭受捕捞型过度捕捞[1]；SS3 和VPA 模型认为其处于上述两种过度捕捞状态[4,12]；而ASPM 模型则表明其均不处于过度捕捞状态[2]。综合上述模型，ICCAT 认为种群处于资源型过度捕捞而未遭受捕捞型过度捕捞的状态[12]。本研究与之产生差异的原因可能是由状态空间建模与上述几种模型的结构差异及使用的先验信息不同所导致的。本研究得到大西洋黄鳍金枪鱼的环境容纳量K为178×104t，内禀增长率r为0.210，与同样基于剩余产量理论构建的ASPIC 模型的结果相似[1]，说明评估结果较为可信；与ASPIC 模型相比，评估的种群状态脱离了资源型过度捕捞，这可能是因为近两年捕捞死亡率的下降，使黄鳍金枪鱼资源有机会得到部分恢复。

图2 大西洋黄鳍金枪鱼JABBA 模型S1−S8 方案的CPUE 指数趋势Fig.2 Time-series of input CPUE of Atlantic yellowfin tuna and predicted CPUE of S1−S8 scenarios in JABBA

表4 大西洋黄鳍金枪鱼JABBA 模型S1−S8 方案的拟合效果Table 4 Goodness of fitting of S1−S8 scenarios in JABBA for Atlantic yellowfin tuna

图3 大西洋黄鳍金枪鱼JABBA 基础模型参数先验分布（深色）和后验分布（浅色）Fig.3 Priors (dark) and posteriors (light) of parameters of base case in JABBA for Atlantic yellowfin tuna

相对ICCAT 当前的资源评估模型而言，JABBA作为剩余产量模型的一种，其结果可靠性较高但无法充分利用鱼类的生物学数据，与ICCAT 使用的其他剩余产量模型如ASPIC 等相比，JABBA 可以估计过程误差，对形状参数的估计更自由，但JABBA 的基本假设为渔获量不存在误差，这一点有待改进。此外，在v1.1 版本的JABBA 中，在选择Pella-Tomlinson 产量函数时，模型无法将m作为未知参数直接估算，必须通过假定BMSY与K的关系得到，因此下一步我们将探寻JABBA 的更高版本，以研究m的估算问题。JABBA 模型参数设定自由，拟合快速，当前已有用其评估大西洋剑鱼（Xiphias gladius）[23]、大西洋大眼金枪鱼（Thunnus obesus）[24]、大西洋蓝枪鱼（Makaira nigric-ans）[25]等的研究，相信其在RFMOs 的资源评估中将发挥越来越重要的作用。

表5 大西洋黄鳍金枪鱼JABBA 基础模型参数后验估计值及其95%置信区间Table 5 Posterior estimates and 95% confidence intervals of parameter of base case in JABBA for Atlantic yellowfin tuna

近年来大西洋黄鳍金枪鱼的渔获量为13×104t 左右，预测分析结果表明，当前TAC（11×104t）管理措施对其种群的养护是有效的，可以实现ICCAT 的管理目标，而在当前的渔获量水平下，种群生物量仍能保持一定的速度增长。近年来，我国的渔获量仅为0.05×104t 左右[26]，占总渔获量比例较小，且都来自于延绳钓渔业的兼捕渔获，因此我国的大西洋黄鳍金枪鱼渔业仍有一定的开发空间。

大西洋黄鳍金枪鱼渔业主要有围网、延绳钓和饵钓3 种，其中围网渔业占总渔获量的70%左右，且围网渔业主要在东部大西洋作业[26]。而黄鳍金枪鱼的生长分为两个阶段[27-28]，幼鱼阶段生长较为缓慢，成鱼阶段生长快速，且黄鳍金枪鱼幼鱼主要在大西洋东部完成早期生活史[29]。当前围网渔业渔获量上升[26]，造成的黄鳍金枪鱼幼鱼死亡率偏高[12]，可能导致补充量不足，种群内禀增长率降低，致使剩余产量减少，可以考虑适当限制围网渔业的捕捞投入，以更好地养护大西洋黄鳍金枪鱼资源。

图4 1950−2017 年大西洋黄鳍金枪鱼JABBA 模型基础模型资源开发状态变化图Fig.4 Kobe phase plot showing estimated trajectories(1950−2017) of B/BMSYand F/FMSYof Atlantic yellowfin tuna of base case in JABBA

图5 大西洋黄鳍金枪鱼JABBA 基础模型1950−2017 年F/FMSY和B/BMSY趋势Fig.5 F/FMSYand B/BMSYof Atlantic yellowfin tuna from 1950 to 2017 of base case in JABBA

图6 不同TAC 目标下的大西洋黄鳍金枪鱼JABBA 模型基础模型B/BMSY预测（2019−2027 年）Fig.6 Future projection (2019−2027) of B/BMSYof Atlantic yellowfin tuna of base case in JABBA under different TACs

剩余产量模型将种群所有个体生命史的动态变化过程进行了高度综合，模型具有参数少、所需数据相对简单的特点，而形状参数较难准确估计且容易导致资源评估的失败[30]，因此在模型拟合结果相差不大时，本研究最终选择了较简单的Fox 而放弃形状参数不易估计的Pella-Tomlinson 产量函数形式。贝叶斯方法把经验判断、前人的研究结果与现有数据相结合[8,31]，后验概率分布由先验概率分布和模型数据共同决定，但如果所用的数据不包含足够的信息，那么后验概率分布可能完全由先验概率主导和控制[32]。因此在使用贝叶斯资源评估方法时，对后验概率分布与先验概率分布进行比较分析显得尤为重要[33-34]。本研究中的敏感性分析显示，K、r的后验分布对先验分布是否有信息并不敏感，说明数据为模型的贝叶斯方法提供了足够的信息。

回溯性误差在渔业资源评估中比较普遍，误差过大可能导致渔业管理的失败[35]。对基础模型的回溯性分析中，对生物量的估计过低而对捕捞死亡率估计过高，这可能是由近年来黄鳍金枪鱼渔获量下降导致。4 个参数的ρ值均趋近于0，结合图形绘制结果可以表明不存在明显的回溯性误差，这可能是由于状态空间建模不仅给出了传统模型的点估计值，同时能量化观测误差和过程误差的不确定性，从而避免了一定的回溯性问题[20,36]。

本研究表明，早期渔获量数据误报率会对资源量和捕捞死亡率的结果产生一定影响，而种群状态并没有明显改变，不会影响对种群健康状态的判断。一般来说渔获量数据以少报居多，本研究表明此时资源评估的结果将更加乐观。但本研究未考虑其他时间段内渔获量数据失真问题，而近期渔获量数据对当前资源状态的判断有更大影响，此外，下一步的研究还应考虑近期数据的随机误差等情况[37]。

表6 不同TAC 目标下大西洋黄鳍金枪鱼2019−2027 年B>BMSY的概率Table 6 The probability that B>BMSYof Atlantic yellowfin tuna under different TAC targets in 2019−2027

表7 不同TAC 目标下大西洋黄鳍金枪鱼2019−2027 年F>FMSY的概率Table 7 The probability that F>FMSYof Atlantic yellowfin tuna under different TAC targets in 2019−2027

表8 不同TAC 目标下大西洋黄鳍金枪鱼2019−2027 年处于健康状态的概率Table 8 The probability that the Atlantic yellowfin tuna is in healthy status under different TAC targets in 2019−2027

表9 不同渔获量误报比例下大西洋黄鳍金枪鱼JABBA 基础模型评估资源状态Table 9 Stock status of Atlantic yellowfin tuna in different mis-reported rates of catches of base case in JABBA

图7 大西洋黄鳍金枪鱼JABBA 基础模型B、B/BMSY、F、F/FMSY的回溯性分析Fig.7 Retrospective analysis of B,B/BMSY,F,F/FMSYof base case in JABBA of Atlantic yellowfin tuna

致谢：感谢渔业资源和生态系统量化评估与管理研究室赵蓬蓬等师兄在论文修改方面的帮助。