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

在JABBA软件中，需要估计的参数包括：承载能力(K)、内禀增长率(r)、初始资源量与K的比值(φ)、形状参数(m)，日本延绳钓渔业的捕捞系数(q_J)、中国台湾延绳钓渔业的捕捞系数(q_J)、处理误差方差(σ_η²)，日本标准化CPUE观察误差方差(σ_{est, J}²)、中国台湾标准化CPUE观察误差方差(σ_{est, T}²)。JABBA软件的具体描述可参见文献[11]。

1.2.2 模型参数的先验设置

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

(2) σ_η²，σ_{est, J}²与σ_{est, T}²的先验设置    假设σ_η²，σ_{est, J}²与σ_{est, T}²的倒数即1/σ_η²、1/σ_{est, J}²与1/σ_{est, T}²服从无信息的伽马分布, 即其形状参数为0.001, 尺度参数的倒数为0.001，记作GM(0.001, 0.001)。

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

1.2.3 其它参数的设置

由于m较难估计，为此，本文假设m分别取0.75、1.0 (即Fox模型)^[12]、1.5、2.0(即Schaefer模型)^[12]及2.5。此外，本文将标准化CPUE固定方差(σ_fix²)设为0.01，两个标准化CPUE外部估计方差设为0。

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

1.2.4 模型收敛诊断

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

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

$ DIC = 2\bar D(s, \theta ) - D(s, \bar \theta )。$

(1)

式中：θ为参数；s为数据；θ为参数θ的后验平均值；D(s, θ)为后验平均偏差(Posterior mean deviance)；D(s, θ)为后验均值偏差，可由式(2)计算。

$ D(s, \theta ) = - 2{\rm{ln}}({\rm{L}}(s|\theta ))。$

(2)

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

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

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

$ {P_y} = \left( {{P_{y - 1}} + \frac{r}{{m - 1}}{P_{y - 1}}\left( {1 - {P_{y - 1}}^{m - 1}} \right) - \frac{{{C_{y - 1}}}}{K}} \right){{\rm{e}}^{\eta y}}; $

(3)

$ {P_y} = \frac{{{B_y}}}{K}; $

(4)

$ {{\hat I}_{i, y}} = {q_i}{B_y}。$

(5)

式中：i可以为T或J，分别表示中国台湾或日本预测CPUE；y为年份；η_y为正态分布随机数，该正态分布的方差为σ_η²，均值为0；${{\hat I}_{i, y}}$为预测的CPUE。

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

1.3.2 预测能力评价

利用2006—2014年标准化CPUE(I_{i, y, O})及模型预测CPUE(Î_{i, y, P})数据，计算去中心均方根误差E_i：

$ {E_i} = \sqrt {\frac{1}{9}\sum\limits_{y = 2006}^{2014} {{{\left[ {\left( {{I_{i, y, O}} - \overline {{I_{i, O}}} } \right) - \left( {{{\hat I}_{i, y, P}} - \overline {{{\hat I}_{i, P}}} } \right)} \right]}^2}} } ; $

(6)

$ {\overline {{I_{i, O}}} = \frac{1}{9}\sum\limits_{y = 2006}^{2014} {{I_{i, y, O}}} ;} $

(7)

$ {\overline {{I_{i, O}}} = \frac{1}{9}\sum\limits_{y = 2006}^{2014} {{{\hat I}_{i, y, P}}} 。} $

(8)

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

为分析E_i的构成，E_i可由日本或中国台湾标准化CPUE方差(σ_{i, O}²)、模型预测CPUE方差(σ_{i, P}²)及其之间的相关系数(ρ)计算^[15]：

$ {E_i} = \sqrt {\sigma _{i, O}^2 + \sigma _{i, P}^2 - 2{\sigma _{i, O}}{\sigma _{i, P}}{\rho _i}} 。$

(9)

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

$ {{E'}_i} = \sqrt {1 + \delta _{i, P}^2 - 2{\delta _{i, P}}{\rho _i}} ; $

(10)

$ {{E'}_i} = \frac{{{E_i}}}{{{\sigma _{i, O}}}}; $

(11)

$ {\delta _{i, P}} = \frac{{{\sigma _{i, P}}}}{{{\sigma _{i, O}}}}。$

(12)

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

1.4 管理风险分析

在渔业管理中，通常根据捕捞控制规则及评估模型投影的结果确定最终TAC。为评估管理风险，本文利用不同时段(即1950—2005年与1950—2014年)数据估计的参数及式(3)与式(4)，分别以2006年为起始年，进行投影计算TAC，以比较不同时段数据、不同模型估计的TAC的差异。确定TAC的捕捞控制规则是：在2014年，当年资源量(B₂₀₁₄)大于最大可持续产量下的资源量(B_MSY)、并且当年捕捞死亡系数(F₂₀₁₄)小于最大可持续产量下的捕捞死亡系数(F_MSY)的概率大于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拟合较好。

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

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

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)。

2.3 模型估计的TAC比较

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

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的值可由最大可持续产量下的产卵生物量(SSB_MSY)与初始条件下的产卵生物量(SSB₀)的比值确定。由于SSB_MSY受选择系数等具有时变特点的参数影响，从这个角度看，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必须使捕捞死亡系数小于F_MSY，资源量(或产卵生物量)大于B_MSY的概率大于一个确定数，这个数通常是50%或60%。就本文而言，当使用1950—2005年数据时，DIC选择的最佳模型S6的TAC分别为4.77×10⁵ t (50%)和4.58×10⁵ t (60%)，而使用1950—2014年数据时，最佳模型S9的TAC分别为3.22×10⁵ t (50%)和3.18×10⁵ t (60%)，这两者估计的TAC相差1.4×10⁵ t以上，这将使渔业管理存在极大风险。因此，通过后向预报方法不仅可发现DIC选择的最佳模型不能提供最佳预测、DIC选择的最佳模型不能保持时间上的一致性，并可揭示这些模型提供的TAC存在巨大差异，这将使渔业管理者更加了解渔业资源评估的不确定性及其可能引起的渔业管理风险。当需要利用这些剩余产量模型为印度洋黄鳍金枪鱼的渔业管理提供TAC时，后向预报方法提供的结果将提醒渔业管理者必须考虑这些不确定性。如本文，渔业管理者应将F₂₀₁₄ < F_MSY且B₂₀₁₄>B_MSY的概率提高，至少应大于60%以上，以增加渔业管理的稳健性、确保TAC能使渔业管理目标得以实现。