上海海洋大学学报  2019, Vol. 28 Issue (2): 298-304    PDF    
剩余产量模型形状参数对印度洋黄鳍金枪鱼资源评估的影响
官文江1,2, 吴佳文1     
1. 上海海洋大学 海洋科学学院, 上海 201306;
2. 大洋渔业资源可持续开发省部共建教育部重点实验室, 上海 201306
摘要:剩余产量模型中的Schaefer模型与Fox模型在渔业资源评估中被广泛应用,但这两个模型是剩余产量模型一般形式Pella-Tomlinson模型的两个特例,分别由形状参数的两个不同值确定。由于形状参数直接影响种群的生产能力,并与种群的年龄结构、繁殖能力等紧密相关,将形状参数固定为某个值,可能会影响剩余产量模型评估结果的可靠性,为此,本文利用印度洋黄鳍金枪鱼(Thunnus albacares)数据,分析剩余产量模型形状参数对渔业资源评估的影响。结果表明:(1)形状参数较难估计,并随数据时段的不同,形状参数的最佳取值范围会有很大的不同;(2)形状参数会对承载能力、内禀增长率的估计产生显著影响;(3)随形状参数值的增加,资源被掏空率与过度捕捞程度会不断增加,因此,不同形状参数的设置会对渔业资源状态及过度捕捞程度的判断产生重要影响。
关键词剩余产量模型    形状参数    资源评估    黄鳍金枪鱼    印度洋    

尽管剩余产量模型因忽略种群的年龄或体长结构等信息而受到批评[1],但其具有参数与数据需求少、容易理解与执行、直观等特点,同时能提供渔业管理所需的最大可持续产量(maximum sustainable yield, MSY)等生物参考点[2],在渔业资源评估中一直被广泛应用,如在大西洋大眼金枪鱼(Thunnus obesus)[3]、印度洋黄鳍金枪鱼(Thunnus albacares)[4]、印度洋长鳍金枪鱼(Thunnus alalunga)[5]、南大西洋剑鱼(Xiphias gladius)[2]等资源评估中被广泛应用,在最近开发的渔业管理策略评价(management strategy evaluation, MSE)中,也因为其具有简单等特点而被广泛用作操作模型或评估模型[6]

PELLA与TOMLINSON[7]给出了剩余产量模型的一般形式,即Pella-Tomlinson模型[8]。在Pella-Tomlinson模型下,形状参数取值不同则可获得不同的剩余产量模型。其中,以Schaefer模型与Fox模型[8-9]最为有名,由于形状参数较难估计,因而这两个模型被广泛应用。由于Schaefer模型的产量函数具有对称性,该特点常被认为与种群动态实际不符[10],因此,Fox模型的评估结果更容易被接受[5, 10]。形状参数直接影响种群的生产能力,并与种群的年龄结构、繁殖能力等紧密相关,将形状参数固定为某个值(如Fox模型),则会影响剩余产量模型评估结果的可靠性。目前,对形状参数的设置及其所带来的影响的讨论较少。为此,本文尝试利用印度洋黄鳍金枪鱼数据,分析剩余产量模型形状参数对渔业资源评估的影响,以期为科学使用剩余产量模型、避免渔业管理风险提供科学依据。

1 材料与方法 1.1 数据来源

数据来自印度洋金枪鱼委员会(Indian Ocean Tuna Commission, IOTC)(http://www.iotc.org/meetings/17th-working-party-tropical-tunas-wptt17),包括渔获量与标准化CPUE(catch per unit effort)数据。渔获量数据的时间跨度为1950—2014年。标准化CPUE数据分别来自中国台湾与日本延绳钓渔业。中国台湾标准化CPUE数据的时间跨度为1980—2012年,而日本标准化CPUE数据的时间跨度为1963—2014年。由于1971年及以前的日本标准化CPUE数据不能真实反映黄鳍金枪鱼资源量的变化,因此该时段的日本标准化CPUE数据被舍弃[4, 11]

1.2 剩余产量模型

Pella-Tomlinson模型的剩余产量函数[2]如下:

    (1)
    (2)

式中:m为形状参数,Byy年的生物量,Syy年的剩余产量,K为承载能力,r为内禀增长率, Cyy年的总渔获量。当m为2时,式(2)为Schaefer模型;当m趋于1时,式(2)为Fox模型[8-9]

Pella-Tomlinson模型的MSY及在MSY下的生物量(BMSY)、捕捞死亡系数(FMSY)与rmK有如下关系:

    (3)
    (4)
    (5)
1.3 参数估计 1.3.1 参数估计软件

采用JABBA(just another bayesian biomass assessment)软件[2]进行参数估计,该软件的种群动态方程在式(2)的基础上增加了处理误差,并重新参数化为式(6):

    (6)
    (7)

式中:ηy为处理误差,其服从方差为ση2,均值为0的正态分布,当y=1时,φ=P1

JABBA软件的观测模型为

    (8)

式中:i可以为T或J,分别表示中国台湾或日本标准化CPUE,q为捕捞系数。εy, i为观测误差,其服从方差为σε, y, i2,均值为0的正态分布。σε, y, i2由3部分组成:

    (9)

式中:为外部估计方差(如CPUE标准化时估计的方差),本文设为0;σfix2为固定方差,本文设为0.01,σest, i2由JABBA估计。因此,JABBA需要估计的参数集为θ={K, r, φ, m, qT, qJ, ση2, σest, J2, σest, T2}。JABBA软件的具体描述可参见文献[2]。

1.3.2 先验设置

(1) rKφqJqT的先验设置

依官文江等[4]r服从中值与变异系数分别为0.46与0.22的对数正态分布时,模型能获得较好的估计效果,因此,本文采用该先验分布。根据多年资源评估结果[4],为足以覆盖K的合理范围,本文假设K的最小值为最大渔获量(5.29×105 t)的2倍(即1.06×106 t),最大值为最大渔获量的32倍(即1.69×107 t),因此,K服从1.06×106 t至1.69×107 t的均匀分布。由于早期渔获量较低,并且标准化CPUE数据始于1972年,数据无估计φ的有效信息、φ在0至1之间的值均能满足模型[4],为此,本文假设φ服从中值为0.90、变异系数为0.001的对数正态分布。基于标准化CPUE与资源量的正比假设,标准化CPUE对应的捕捞系数一般小于1.0,为足以覆盖捕捞系数值的范围,本文采用JABBA的默认假设,即qJqT均服从0至100的均匀分布。

(2) ση2σest, J2σest, T2的先验设置

假设ση2σest, J2σest, T2的倒数即1/ση2、1/σest, J2与1/σest, T2服从无信息的伽马分布, 即其形状参数为0.001, 尺度参数的倒数为0.001。

1.3.3 计算场景、模型收敛诊断及模型选择

将渔获量与标准化CPUE数据按不同年份分为4个数据集,即1950—2014年、1950—2011年、1950—2008年和1950—2005年等4个数据集。

对上述每个数据集,依次将m假设为0.102 564 1、0.205 128 2、0.307 692 3、…、3.999 9999(即以0.102 564 1为等差的数列,共39个点)等,并利用JABBA估计其他参数的中位数,以分析:1) m可能的最佳取值范围;2) m的取值对资源状态判断的影响; 3) m的最佳取值范围是否会随数据时段的不同而不同。

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

对每个数据集,由m的不同,共可产生39个模型。为确定m的最佳取值范围,本文采用偏差信息准则[13](deviance information criterion, DIC)确定m的最佳取值,即DIC越小,模型拟合越好,该模型对应的m取值也最佳。为考虑模型拟合与模型选择存在的不确定性,本文先采用局部多项式回归模型(local polynomial regression model, LPRM)拟合m与DIC的关系,然后利用LPRM预测不同m对应的DIC值,并以最小DIC、最小DIC+4为范围[14]确定m的最佳取值范围。

2 结果

当使用1950—2014年或1950—2011年数据,且m在1.37至1.92或1.50至2.10时,模型能获得相对较小的DIC(图 1)。但当使用1950—2008年或1950—2005年数据时,要使模型获得相对较小的DIC,m的取值应小于1.0 (图 1)。因此,在使用不同时段数据时,为使模型取得相对较小的DIC,m的取值范围存在较大差异。

图 1 形状参数与DIC的关系 Fig. 1 Relationship between the shape parameter and DIC

图 2可知,当使用1950—2014年数据(其他数据集类似)时,BMSYK的比值随m值的增大而增大,但r的估计随m值的增大则是先增大,并在m为1.95时,r取得最大值,此后随m值的增加,r的值逐渐减少;K的估计值随m值的增加而单调增加;MSY随m的增加而具有减少趋势,但在m小于1.0时,单调下降幅度相对较大,在m大于1.0、小于1.5时,MSY随m值的增加几乎不变,而在m大于1.5时,MSY又随m值的增加开始缓慢下降。

图 2 形状参数对BMSYK比值、rK及MSY的影响 Fig. 2 The impacts of shape parameter on the ratio of BMSY to K, r, K and MSY

图 3可知, 当使用1950—2014年数据时,m小于0.41时,最后1年(2014年)的生物量(即B2014)与BMSY的比值随m的增加而大幅减少,但当m大于0.41时,该比值随m的增加而减少的幅度变得相对平缓;最后1年的捕捞死亡系数(即F2014)与FMSY的比值随m的增加而增加,当m小于0.72时,F2014FMSY的比值增加幅度较大,但当m大于0.72时,该比值的增加幅度变得相对平缓;除m小于0.41的情形外,m的变化对F2014FMSY的比值影响较大,而对B2014BMSY的比值影响较小。当使用1950—2005数据(其他数据时段类似,为简化,在此省略)时,也存在类似特点,不同点在于F2005FMSY的比值随m的增加,其增幅更为平缓(图 3)。

图 3 形状参数对资源状态估计的影响 Fig. 3 The impacts of shape parameter on the estimates of the stock status
3 讨论 3.1 标准化CPUE的选择及先验、权重的设置

利用贝叶斯剩余产量模型对印度洋黄鳍金枪鱼进行资源评估,需讨论不同标准化CPUE、先验及数据权重设置等影响,以选择最合适的模型配置进行资源评估。但本文仅利用了官文江等[4]的最佳模型配置设置了模型,以着重讨论形状参数的影响,对标准化CPUE选择及先验、权重设置的相关讨论可参考相关文献。

3.2 形状参数的估计

MCALLISTER等[15]重新参数化了Pella-Tomlinson模型,提出Fletcher-Schaefer模型,并利用mr与种群统计参数(demographic parameter)的复杂经验关系构建了mr的联合先验分布,以提高资源评估质量。从其结果可知,m的后验分布与先验分布基本一致,因此,资源评估数据没有包含估计m的足够信息[15]。LI等[16]在印度洋长鳍金枪鱼的资源评估中也估计了m,但其估计的初始年资源量与次1年资源量相差较大的结果表明该模型不一定合理。

本文结果也表明m较难被估计。如当使用1950—2014年数据时,尽管平滑函数能较好地给出m的估计范围(图 1),但对每个具体的m,其对应的DIC仍存在较大波动,使其全局极小值较难准确获取,如m在1.37至1.92区间能使DIC小于-298.0(即最小DIC加4的值),但在m大于1.92时(如2.05、2.15与2.56)也能使DIC小于-298.0(图 1)。当使用1950—2008或1950—2005年数据时,m对应最小DIC的区间则存在较大的变化(图 1),从而无法确定哪个估计更接近真实情况。

由式(3)可知,如果能确定BMSYK的比值,则m就可被确定(图 2)。为此,WINKER等[2]假设年龄结构模型估计的最大可持续产量下的产卵生物量(SSBMSY)与初始条件下的产卵生物量(SSB0)的比值等于BMSYK的比值,以确定m的值。若采用该假设,多年印度洋黄鳍金枪鱼资源评估结果[11, 17-20]表明,m的取值在0.20~1.0均有分布,但无明显时间趋势。该结果与利用1950—2005年或1950—2008年数据集估计的结果相似(图 1),而与利用1950—2011年或1950—2014年数据集估计的结果不同(图 1)。由于年龄结构模型估计的SSBMSY与SSB0的比值不一定可靠,同时该比值和BMSYK比值的确切关系仍有待进一步研究,因此无法推断该结果是否正确。但这一假设为获取m提供了新思路,更为重要的是,由于SSBMSY与选择系数紧密关联,并随选择系数的改变而改变,因此,可以推定m也应该是变化的,这使m的估计更难。

尽管剩余产量模型在渔业资源评估中使用较广,但估计m的实例不多,以致Pella-Tomlinson模型较少被使用[15]。当前估计m可能取值范围的几种常用方法是:(1)利用mr及其与种群统计参数的复杂经验关系确定m;(2)利用m的DIC剖面(或其他模型选择量剖面,如似然值剖面)选择合适的m,以防止获得局部最优估计;(3)利用年龄结构模型估计的SSBMSY与SSB0的比值确定m,而直接对m进行估计则有可能无法获得全局最优估计。

3.3 形状参数对资源评估的影响

图 2图 3表明,m直接影响参数rK的估计,进而影响MSY、最后1年生物量(Blast)与BMSY比值(对1950—2014年的数据,为B2014/BMSY)和最后1年捕捞死亡系数(Flast)与FMSY比值(对1950—2014年的数据,为F2014/FMSY)的估计。而MSY,特别是Blast/BMSYFlast/FMSY是非常重要的管理参数。Blast/BMSY小于1则发生资源型过度捕捞,而Flast/FMSY大于1则出现捕捞型过度捕捞,这两种过度捕捞都可能引起管理措施的调整及总可捕量(total allowable catch, TAC)的下调。并且,由图 3可知:随m的增大,资源评估结果得出资源掏空程度与遭受的捕捞压力是不断增加的;当m小于0.41时,m的增加对资源掏空程度的判断影响较大;当m大于0.41时,m的增加则对捕捞型过度捕捞程度的判断影响较大。因此,m设置不当将影响对资源状态及过度捕捞程度的判断。在渔业管理中,应采用渔业管理策略评价等方法以避免这些不确定性的影响[21-22]

此外,由图 2可知,随m的不同,r的估计不同。但若r为内禀增长率,依据其与种群统计参数的关系[21]r不应随m而改变。因此,当使用Pella-Tomlinson模型时,除m取2之外,参数r不再具有内禀增长率的内在含义[15],在r的先验设置时应注意这一点或采用其他参数化模型,如MCALLISTER等[15]提出的Fletcher-Schaefer模型,这类模型能使r一直满足内禀增长率的原始定义[8]

参考文献
[1]
WANG S P, MAUNDER M N, AIRES-DA-SILVA A. Selectivity's distortion of the production function and its influence on management advice from surplus production models[J]. Fisheries Research, 2014, 158: 181-193. DOI:10.1016/j.fishres.2014.01.017
[2]
WINKER H, CARVALHO F, KAPUR M. JABBA:just another Bayesian biomass assessment[J]. Fisheries Research, 2018, 204: 275-288. DOI:10.1016/j.fishres.2018.03.010
[3]
KELL L, MERINO G. Stock assessment diagnostics for Atlantic Bigeye tuna[J]. ICCAT Collective Volume of Scientific Papers, 2016, 72(1): 245-265.
[4]
官文江, 朱江峰, 田思泉. 应用贝叶斯生物量动态产量模型评估印度洋黄鳍金枪鱼资源[J]. 中国水产科学, 2018, 25(3): 621-631.
GUAN W J, ZHU J F, TIAN S Q. Assessment of the Indian Ocean yellowfin tuna (Thunnus albacares) by using a Bayesian biomass dynamic model[J]. Journal of Fisheries Sciences of China, 2018, 25(3): 621-631.
[5]
GUAN W J, TANG L, ZHU J F, et al. Application of a Bayesian method to data-poor stock assessment by using Indian Ocean albacore (Thunnus alalunga) stock assessment as an example[J]. Acta Oceanologica Sinica, 2016, 35(2): 117-125. DOI:10.1007/s13131-016-0814-0
[6]
MERINO G, ARRIZABALAGA H, SANTIAGO J, et al. Updated evaluation of harvest control rules for North Atlantic albacore through management strategy evaluation[J]. ICCAT Collective Volume of Scientific Papers, 2017, 74(2): 457-478.
[7]
PELLA J J, TOMLINSON P K. A generalized stock production model[J]. Inter-American Tropical Tuna Commission Bulletin, 1969, 13: 421-458.
[8]
QUINN T J, DERISO R B. Quantitative fish dynamic[M]. New York: Oxford University Press, 1999: 398-401.
[9]
詹秉义. 渔业资源评估[M]. 北京: 中国农业出版社, 1995.
ZHAN B Y. Fishery stock assessment[M]. Beijing: China Agricultural Press, 1995.
[10]
IOTC. 2014. Report of the fifth session of the IOTC working party on temperate tunas. IOTC-2014-WPTmT05-R[E] [R]. IOTC 5th Working Party on Temperate Tunas, Indian Ocean Tuna Commission, Busan, Rep. of Korea, 28-31 July 2014.
[11]
LANGLEY A. Stock assessment of yellowfin tuna in the Indian Ocean using Stock Synthesis (IOTC-2015-WPTT17-30)[R]. 17th Session of the IOTC Working Party on Tropical Tunas, Montpellier, France: IOTC, 22-28 October 2015.
[12]
KÉRY M. Introduction to WinBUGS for ecologists:a Bayesian approach to regression, ANOVA, mixed models and related analyses[M]. San Diego: Academic Press, 2010.
[13]
SPIEGELHALTER D J, BEST N G, CARLIN B P, et al. Bayesian measures of model complexity and fit (with discussion)[J]. Journal of the Royal Statistical Society:-Series B (Statistical Methodology), 2002, 64(4): 583-639. DOI:10.1111/rssb.2002.64.issue-4
[14]
官文江, 田思泉, 王学昉, 等. CPUE标准化方法与模型选择的回顾与展望[J]. 中国水产科学, 2014, 21(4): 852-862.
GUAN W J, TIAN S Q, WANG X F, et al. A review of methods and model selection for standardizing CPUE[J]. Journal of Fisheryies Sciences of China, 2014, 21(4): 852-862.
[15]
MCALLISTER M K, BABCOCK E A, PIKITCH E K, et al. Application of a non-equilibrium generalized production model to South and North Atlantic swordfish:combining Bayesian and demographic methods for parameter estimation[J]. ICCAT Collective Volume of Scientific Papers, 2000, 51: 1523-15501.
[16]
LI B, CAO J, ZHU J. Analyzing population dynamics of Indian Ocean albacore (Thunnus alalunga) using Bayesian state-space production model (IOTC-2016-WPTmT06-24)[R]. Sixth Session of the IOTC Working Party on Temperature Tunas., Shanghai, China: IOTC, 18-21 July 2016.
[17]
SHONO H, SATOH K, OKAMOTO H. Preliminary stock assessment for yellowfin tuna in the Indian Ocean using stock synthesis Ⅱ (SS2) (IOTC-2007-WPTT-11)[R]. 9th Session of the IOTC Working Party on Tropical Tunas, Victoria, Seychelles, 16-20 July 2007.
[18]
LANGLEY A, HERRERA M, MILLION J. Stock assessment of yellowfin tuna in the Indian Ocean using Multifan-CL (IOTC-2010-WPTT-23)[R]. 12th Session of the IOTC Working Party on Tropical Tunas, Victoria, Seychelles, 18-25 October 2010.
[19]
LANGLEY A, HERRERA M, MILLION J. Stock assessment of yellowfin tuna in the Indian Ocean using MULTIFAN-CL (IOTC-2012-WPTT14-38)[R]. IOTC, 14th Session of the IOTC Working Party on Tropical Tunas, Mauritius, 24-29 October 2012.
[20]
LANGLEY A. An update of the 2015 Indian Ocean yellowfin tuna stock assessment for 2016(IOTC-2016-WPTT18-27)[R]. 18th Session of the IOTC Working Party on Tropical Tunas, Seychelles, 5-10 November 2016.
[21]
李志, 王家启, 田思泉. 印度洋黄鳍金枪鱼渔业管理策略评价的初步研究[J]. 上海海洋大学学报, 2016, 25(2): 255-262.
LI Z, WANG J Q, TIAN S Q. Management strategy evaluation for yellowfin tuna fishery in Indian Ocean[J]. Journal of Shanghai Ocean University, 2016, 25(2): 255-262.
[22]
官文江, 高峰, 陈新军. 卫星遥感在海洋渔业资源开发、管理与保护中的应用[J]. 上海海洋大学学报, 2017, 26(3): 440-449.
GUAN W J, GAO F, CHEN X J. Review of the application of satellite remote sensing in the exploitation, management and protection of marine fisheries resources[J]. Journal of Shanghai Ocean University, 2017, 26(3): 440-449.
[23]
MCALLISTER M K, PIKITCH E K, BABCOCK E A. Using demographic methods to construct Bayesian priors for the intrinsic rate of increase in the Schaefer model and implications for stock rebuilding[J]. Canadian Journal of Fisheries and Aquatic Sciences, 2001, 58(9): 1871-1890. DOI:10.1139/f01-114
Impacts of shape parameter of surplus production model on stock assessment of Indian Ocean yellowfin tuna
GUAN Wenjiang1,2, WU Jiawen1     
1. College of Marine Sciences, Shanghai Ocean University, Shanghai 201306, China;
2. The Key Laboratory of Sustainable Exploitation of Oceanic Fisheries Resources, Ministry of Education, Shanghai 201306, China
Abstract: Schaefer form and Fox form of surplus production model are widely used in fishery stock assessment. However, the two models are two special cases of the generalized surplus production model, i.e. Pella-Tomlinson model, with two special values of the shape parameter. Because the shape parameter has a close relationship with the age structure and reproductive capacity of the population and directly affects the modeled population productivity, fixing the shape parameter at some special values would raise doubts about the results of the stock assessment models of the model. The impacts of shape parameter on stock assessment results were analyzed based on Indian Ocean yellowfin tuna (Thunnus albacares) data. The results showed:1) the shape parameter was difficult to be estimated and its best value range would be greatly different with different data sets from different time periods; 2) the shape parameter had obvious impacts on the estimate of carrying capacity and the intrinsic rate of population increase; 3) the depletion levels of the stock and the overfished or overfishing levels were enhanced with the value of the shape parameter increasing and accordingly the shape parameter had important impacts on the inference of stock status and overfished and/or overfishing levels.
Key words: Surplus production model     shape parameter     stock assessment     Thunnus albacares     Indian Ocean