利用后向预报方法分析印度洋黄鳍金枪鱼资源评估模型❋

官文江, 吴佳文, 曹友华

(1. 上海海洋大学海洋科学学院，上海 201306； 2. 上海海洋大学大洋渔业资源可持续开发省部共建教育部重点实验室， 上海 201306)

印度洋黄鳍金枪鱼(Thunnusalbacares)是商业金枪鱼渔业的重要捕捞目标鱼种之一，主要被法国、西班牙的围网渔业，印度尼西亚、中国台湾的冰鲜延绳钓渔业，中国台湾、日本、韩国的超低温延绳钓渔业及马尔代夫、印度等国的手钓、竿钓、刺网等渔业所利用[1]。近年，随着印度洋海盗活动减少、安全形式好转，印度洋黄鳍金枪鱼渔获量又得以大幅提升，年平均产量超过了4.0×105t。尽管我国捕捞的印度洋黄鳍金枪鱼较少，但在我国的印度洋金枪鱼渔业中，其产量仅次于大眼金枪鱼(Thunnusobesus)与长鳍金枪鱼(Thunnusalalunga)。印度洋黄鳍金枪鱼的资源状况也日益受到我国学者的关注[1-5]。

当前，剩余产量模型(Bayesian Biomass Dynamic Model，BBDM)[1]、年龄结构模型(Statistical Catch at Age，SCAA)[6]及合成模型(Stock Synthesis，SS3)[7-8]被用于印度洋黄鳍金枪鱼的资源评估。尽管SS3模型的评估结果被用于确定印度洋黄鳍金枪鱼的资源状态[9]，但上述各模型的评估结果均存在较大不确定性[1,9]。特别是这些模型的诊断与选择，主要采用的依据是模型对观察数据(主要是资源丰度指数及年龄或体长组成数据)的拟合度。但对渔业管理而言，如总可捕量(Total Allowable Catch，TAC)的确定等，模型的预测能力也非常重要。在当前的渔业资源评估中，相关学者并未足够重视模型的预测能力，很少评价模型的预测能力、并将其作为评价渔业资源评估与管理质量的依据[10]。为此，本文以印度洋黄鳍金枪鱼资源评估为例，利用后向预报方法评价模型的预测能力，以分析各模型渔业资源评估与管理质量，为科学选择渔业资源评估模型、避免渔业管理风险提供科学依据。

1 材料与方法

1.1 数据来源

用于印度洋黄鳍金枪鱼资源评估的数据包括渔获量数据与标准化CPUE(Catch Per Unit Effort)数据。其中，标准化CPUE数据分别来自中国台湾与日本延绳钓渔业。渔获量数据的时间跨度为1950—2014年，中国台湾标准化CPUE数据的时间跨度为1980—2012年，而日本标准化CPUE数据的时间跨度为1963—2014年[1]。由于1972年以前的日本标准化CPUE数据存在问题[1,7]，因此，本文仅使用了1972年及以后的数据。上述数据均来自印度洋金枪鱼委员会(Indian Ocean Tuna Commission，IOTC)网站(http://www.iotc.org/meetings/17th-working-party-tropical-tunas-wptt17)。

1.2 评估模型及配置

1.2.1 评估模型 本文使用贝叶斯剩余产量模型对印度洋黄鳍金枪鱼进行资源评估，该模型由软件JABBA(Just Another Bayesian Biomass Assessment)实现参数估计[11]。

1.2.2 模型参数的先验设置

(1)K，φ，qJ及qT的先验设置 根据文献[1]，假设K服从1.06×106～1.69×107t的均匀分布，记作U(1.06×106, 1.69×107)，φ服从中值为0.90，变异系数为0.001的对数正态分布，记作L(0.9, 0.001)；假设qJ与qT均服从0至100的均匀分布[11]，记作U(0, 100)。

(3)r的先验设置 根据文献[1]，假设r服从中值与变异系数分别为 0.46与0.22及0.75与0.15的对数正态分布，分别记作L(0.46,0.22)、L(0.75,0.15)。

根据m及r先验设置的不同，共有10个模型(见表1)。

1.2.4 模型收敛诊断 模型收敛诊断采用Gelman-Rubin 统计量，并以1.1为阈值，即Gelman-Rubin 统计量小于1.1，则认为模型收敛[13]。本文仅对收敛模型进行分析。

表1 不同模型的先验设置

注：*表示在实际计算中取1.001，因为1.0会使模型分母为0[11]；其他参数的先验分布在各模型中均相同，因此表1未列出。

Note:*denote the 1.001was used in practical calculations, because the 1.0 make the denominator in the model become zero[11]; the priors for the other parameters are same for the 10 models and are not shown in Table 1.

模型拟合效果结合偏差信息准则(Deviance Information Criterion, DIC)[14]与泰勒图[15]进行分析评价，其中DIC由式(1)计算。

(1)

D(s,θ)=-2ln(L(s|θ))。

(2)

式中：ln为自然对数；L为似然函数。

1.3 后向预报及预测能力评价

1.3.1 后向预报 为执行后向预报，本文将1950—2005年期间的数据用于估计参数，2006—2014年的渔获量数据用于后向预报、并预测日本及中国台湾延绳钓渔业CPUE的值，依JABBA软件[11]，其计算方程如下：

(3)

(4)

(5)

对式(5)多次迭代的结果取中值，则得到日本及中国台湾延绳钓渔业CPUE的预测值。

(6)

(7)

(8)

当Ei越接近0时，模型预测能力越好[10]。

(9)

为更好比较不同数据对观测的预测效果，需去除不同观测方差影响，因此，式(9)等式两边同时除以σi,O，得:

(10)

(11)

(12)

由式(10)可知，当模型预测CPUE方差等于标准化CPUE方差(此时为1)，且相关系数为1时，E′i为0，模型预测能力最好。

1.4 管理风险分析

在渔业管理中，通常根据捕捞控制规则及评估模型投影的结果确定最终TAC。为评估管理风险，本文利用不同时段(即1950—2005年与1950—2014年)数据估计的参数及式(3)与式(4)，分别以2006年为起始年，进行投影计算TAC，以比较不同时段数据、不同模型估计的TAC的差异。确定TAC的捕捞控制规则是：在2014年，当年资源量(B2014)大于最大可持续产量下的资源量(BMSY)、并且当年捕捞死亡系数(F2014)小于最大可持续产量下的捕捞死亡系数(FMSY)的概率大于G。G通常设为50%和60%。

2 结果

2.1 模型拟合效果

由图1可知，模型预测的CPUE与日本标准化CPUE具有较高的相关系数(大于0.8)，与日本标准化CPUE所具有的方差(图1为1.0)相比，模型预测CPUE的方差偏小(小于0.5)；模型预测的CPUE与中国台湾标准化CPUE的相关系数较小(小于0.45)，但模型预测CPUE的方差偏大(大于0.9)。模型预测CPUE与日本标准化CPUE的去中心均方根误差(即图1实心点至空心点[1.0, 0] 的距离)小于0.8，而模型预测CPUE与中国台湾标准化CPUE的去中心均方根误差大于1.0，因此所有模型均对日本标准化CPUE拟合较好。

图1 不同模型对CPUE拟合效果的泰勒图

随形状参数值的增加，模型预测CPUE与日本标准化CPUE的相关系数变化较小，但模型预测CPUE的方差不断减少，使其去中心均方根误差逐渐变大。中国台湾标准化CPUE与模型预测CPUE具有类似关系，但其去中心均方根误差随模型预测CPUE方差变小而逐渐变小(见图1)。因此，随m的增加，模型对日本标准化CPUE的拟合效果变差，而对中国台湾标准化CPUE的拟合效果则逐渐变好。

若采用DIC选择模型，则模型S6或S7应为最佳模型(见表2)。

表2 不同模型的DIC估计的去中心均方根误差、预测标准差与相关系数

注：*表示P<0.05。*denoteP<0.05.

2.2 模型预测能力分析

当采用模型S1～S10进行预测时，所有模型预测的CPUE均大于日本标准化CPUE，并且当模型为S1、S2、S6和S7时，大部分日本标准化CPUE数据点不在模型预测CPUE的95%预测区间内(见图2)。模型预测CPUE与日本标准化CPUE均呈负相关关系，甚至是显著的负相关关系(见表2)。模型预测CPUE与中国台湾标准化CPUE关系略好(见表2、图2)，均呈正相关关系，但不显著(见表2)。同时，由表2可知，随m的增加，模型预测能力有逐渐变好的趋势，即去中心均方根误差逐渐变小(见表2)。并且，当r先验设为L(0.46, 0.22)时(即模型S1～S5)，其模型预测效果好于r先验设为L(0.75, 0.15)时(即模型S6～S10)的模型预测效果，模型S5的预测能力最好，但其DIC最大(见表2)。

(点为标准化CPUE，实线为模型预测CPUE，始于2006年的虚线表示95%预测区间。The points was observed CPUE； The solid lines was predicted CPUE； The dash lines began at 2006 denoted 95% prediction interval.)

图2 观测CPUE与模型预测CPUE

Fig.2 Observed CPUE versus predicted CPUE

2.3 模型估计的TAC比较

使用相同模型、不同时段数据估计的TAC存在较大差异，如当使用模型S6及1950—2005年数据时，在F2014BMSY的概率大于50%和60%的情况下、其TAC分别为4.77×105和4.58×105t，而使用相同模型及1950—2014年数据时，该值分别为3.69×105和3.66×105t，TAC分别减少了23%与21%。使用1950—2005年数据时，采用DIC选择的最佳模型为S6或S7，而在使用1950—2014年数据时，采用DIC选择的最佳模型为S9或S10，四模型计算的TAC差异更大(见表3)，这使渔业管理存在巨大风险。

表3 不同模型估计的TAC

注：*模型分别使用了1950至2005年的渔获量数据及1980—2005年、1972—2005年的中国台湾与日本标准化CPUE数据；**模型分别使用了1950至2014年的渔获量数据及1980至2014年、1972至2014年的中国台湾与日本标准化CPUE数据;***表示模型具有最小DIC。

Note：*denote the catch data from 1950 to 2005 and the standardized CPUE from 1980 to 2005 from China Taiwan or from 1972 to 2005 from Japan were used in the model;**denote the catch data from 1950 to 2014 and the standardized CPUE from 1980 to 2014 from China Taiwan or from 1972 to 2014 from Japan were used in the model;***denote the model with minimum DIC.

3 讨论

3.1 标准化CPUE的选择及先验、形状参数的设置

标准化CPUE的选择及参数先验的设置对资源评估结果有重要影响，但不是本文讨论的重点，为此，本文仅参考了官文江等[1]的结果。本文没有考虑r均匀无信息分布的情形，因为使用该先验分布时，模型无法收敛。

在使用剩余产量模型时，由于形状参数较难被估计，因此在渔业资源评估实践中主要使用Fox模型(m为1)或Schaefer模型(m为2)。同时，由于Schaefer模型的产量函数具有对称性，这使得其评估结果较难被接受，因此造成Fox模型更合理的假象，但事实可能并非如此。为此，本文增加了多种形状参数的设置，以展示形状参数的重要影响。本文结果表明，m对参数、生物参考点的估计有重要影响(见表3)，也影响模型的预测能力(见表2、图2)。同时，使用1950—2005年数据与使用1950—2014年数据估计的最佳模型的m不同，可能暗示m不应该作为常数处理。此外，Winker等[11]认为m的值可由最大可持续产量下的产卵生物量(SSBMSY)与初始条件下的产卵生物量(SSB0)的比值确定。由于SSBMSY受选择系数等具有时变特点的参数影响，从这个角度看，m也应该具有时变特点。

3.2 后向预报方法可用于评价渔业资源评估模型的质量

在当前渔业资源评估中，模型诊断与选择主要依赖于模型对观测数据的拟合度[10]，如采用DIC或均方根误差选择模型，但模型对观测数据拟合的好坏，并不一定能决定模型正确与否。如本文，当使用1950—2005年数据时，所有模型均能较好拟合日本标准化CPUE数据，但模型预测的2006—2014年的CPUE均与相应时段日本标准化CPUE呈负相关关系，模型预测能力较差。同样，采用DIC选择的最佳模型S6或S7，其预测能力也最差；并且当使用1950—2014年数据时，采用DIC选择的最佳模型为S9或S10，这表明随着使用数据时段的不同，最佳模型缺少一致性(见表2、表3)。而模型预测能力较差、最佳模型缺少一致性等均表明模型可能没有正确模拟种群的演化动态、渔业资源评估结果可能存在问题。

由于缺少对渔业资源量进行直接观测的有效手段，渔业资源评估与预测结果很难得到验证，而后向预报方法为评价渔业资源评估模型的预测能力与评估质量提供了有效手段，从而有利于渔业资源评估模型的正确诊断与选择。如本文，利用后向预报方法，可评价不同模型的预测能力，评价DIC选择模型的稳定性，从而可在一定程度上判断模型对种群演化动态模拟是否正确，模型评估结果是否存在问题。

3.3 后向预报方法可揭示评估结果的不确定性及其引起的渔业管理风险

在当前金枪鱼渔业管理中，主要采用的捕捞控制规则是：经过若干年 (如10年) 的模型投影后，采用的TAC必须使捕捞死亡系数小于FMSY，资源量(或产卵生物量)大于BMSY的概率大于一个确定数，这个数通常是50%或60%。就本文而言，当使用1950—2005年数据时，DIC选择的最佳模型S6的TAC分别为4.77×105t (50%)和4.58×105t (60%)，而使用1950—2014年数据时，最佳模型S9的TAC分别为3.22×105t (50%)和3.18×105t (60%)，这两者估计的TAC相差1.4×105t以上，这将使渔业管理存在极大风险。因此，通过后向预报方法不仅可发现DIC选择的最佳模型不能提供最佳预测、DIC选择的最佳模型不能保持时间上的一致性，并可揭示这些模型提供的TAC存在巨大差异，这将使渔业管理者更加了解渔业资源评估的不确定性及其可能引起的渔业管理风险。当需要利用这些剩余产量模型为印度洋黄鳍金枪鱼的渔业管理提供TAC时，后向预报方法提供的结果将提醒渔业管理者必须考虑这些不确定性。如本文，渔业管理者应将F2014BMSY的概率提高，至少应大于60%以上，以增加渔业管理的稳健性、确保TAC能使渔业管理目标得以实现。