2. 大洋渔业资源可持续开发省部共建教育部重点实验室, 上海 201306
尽管剩余产量模型因忽略种群的年龄或体长结构等信息而受到批评[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]如下:
式中:m为形状参数,By为y年的生物量,Sy为y年的剩余产量,K为承载能力,r为内禀增长率, Cy为y年的总渔获量。当m为2时,式(2)为Schaefer模型;当m趋于1时,式(2)为Fox模型[8-9]。
Pella-Tomlinson模型的MSY及在MSY下的生物量(BMSY)、捕捞死亡系数(FMSY)与r、m、K有如下关系:
采用JABBA(just another bayesian biomass assessment)软件[2]进行参数估计,该软件的种群动态方程在式(2)的基础上增加了处理误差,并重新参数化为式(6):
式中:ηy为处理误差,其服从方差为ση2,均值为0的正态分布,当y=1时,φ=P1。
JABBA软件的观测模型为
式中:i可以为T或J,分别表示中国台湾或日本标准化CPUE,q为捕捞系数。εy, i为观测误差,其服从方差为σε, y, i2,均值为0的正态分布。σε, y, i2由3部分组成:
式中:
(1) r、K、φ、qJ及qT的先验设置
依官文江等[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的默认假设,即qJ与qT均服从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的取值范围存在较大差异。
由图 2可知,当使用1950—2014年数据(其他数据集类似)时,BMSY与K的比值随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值的增加开始缓慢下降。
由图 3可知, 当使用1950—2014年数据时,m小于0.41时,最后1年(2014年)的生物量(即B2014)与BMSY的比值随m的增加而大幅减少,但当m大于0.41时,该比值随m的增加而减少的幅度变得相对平缓;最后1年的捕捞死亡系数(即F2014)与FMSY的比值随m的增加而增加,当m小于0.72时,F2014与FMSY的比值增加幅度较大,但当m大于0.72时,该比值的增加幅度变得相对平缓;除m小于0.41的情形外,m的变化对F2014与FMSY的比值影响较大,而对B2014与BMSY的比值影响较小。当使用1950—2005数据(其他数据时段类似,为简化,在此省略)时,也存在类似特点,不同点在于F2005与FMSY的比值随m的增加,其增幅更为平缓(图 3)。
利用贝叶斯剩余产量模型对印度洋黄鳍金枪鱼进行资源评估,需讨论不同标准化CPUE、先验及数据权重设置等影响,以选择最合适的模型配置进行资源评估。但本文仅利用了官文江等[4]的最佳模型配置设置了模型,以着重讨论形状参数的影响,对标准化CPUE选择及先验、权重设置的相关讨论可参考相关文献。
3.2 形状参数的估计MCALLISTER等[15]重新参数化了Pella-Tomlinson模型,提出Fletcher-Schaefer模型,并利用m、r与种群统计参数(demographic parameter)的复杂经验关系构建了m与r的联合先验分布,以提高资源评估质量。从其结果可知,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)可知,如果能确定BMSY与K的比值,则m就可被确定(图 2)。为此,WINKER等[2]假设年龄结构模型估计的最大可持续产量下的产卵生物量(SSBMSY)与初始条件下的产卵生物量(SSB0)的比值等于BMSY与K的比值,以确定m的值。若采用该假设,多年印度洋黄鳍金枪鱼资源评估结果[11, 17-20]表明,m的取值在0.20~1.0均有分布,但无明显时间趋势。该结果与利用1950—2005年或1950—2008年数据集估计的结果相似(图 1),而与利用1950—2011年或1950—2014年数据集估计的结果不同(图 1)。由于年龄结构模型估计的SSBMSY与SSB0的比值不一定可靠,同时该比值和BMSY与K比值的确切关系仍有待进一步研究,因此无法推断该结果是否正确。但这一假设为获取m提供了新思路,更为重要的是,由于SSBMSY与选择系数紧密关联,并随选择系数的改变而改变,因此,可以推定m也应该是变化的,这使m的估计更难。
尽管剩余产量模型在渔业资源评估中使用较广,但估计m的实例不多,以致Pella-Tomlinson模型较少被使用[15]。当前估计m可能取值范围的几种常用方法是:(1)利用m、r及其与种群统计参数的复杂经验关系确定m;(2)利用m的DIC剖面(或其他模型选择量剖面,如似然值剖面)选择合适的m,以防止获得局部最优估计;(3)利用年龄结构模型估计的SSBMSY与SSB0的比值确定m,而直接对m进行估计则有可能无法获得全局最优估计。
3.3 形状参数对资源评估的影响图 2与图 3表明,m直接影响参数r与K的估计,进而影响MSY、最后1年生物量(Blast)与BMSY比值(对1950—2014年的数据,为B2014/BMSY)和最后1年捕捞死亡系数(Flast)与FMSY比值(对1950—2014年的数据,为F2014/FMSY)的估计。而MSY,特别是Blast/BMSY,Flast/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 |
2. The Key Laboratory of Sustainable Exploitation of Oceanic Fisheries Resources, Ministry of Education, Shanghai 201306, China