三疣梭子蟹(Portunus trituberculatus)隶属甲壳纲十足目梭子蟹科梭子蟹属，为广温广盐性种类，是中国最重要的海洋经济蟹类之一^[1]。主要分布于日本、朝鲜、菲律宾、马来群岛、红海及中国的渤海、黄海、东海和南海。在东海三疣梭子蟹主要分布在大沙渔场、长江口渔场和舟山渔场^[2]，其渔获量约占全国产量的50%^[3]。东海北部海域是东海最大的三疣梭子蟹产区，其渔获量占东海总渔获量的40%以上^[4]。国外学者对三疣梭子蟹的相关研究较少，Ariyama H^[5]基于1984—2008年日本大阪湾区域三疣梭子蟹的渔获量数据，探究了过于频繁的渔业开发和偶发的人类活动，以及气候事件对三疣梭子蟹补充量造成的影响。我国有关三疣梭子蟹的研究主要集中于养殖领域，包括三疣梭子蟹的人工育苗^[6]、池塘养殖^[7]、常见病害防控^[8]、配合饲料开发^[9-10]以及优良品种选育^[11]等。对野生三疣梭子蟹资源的相关研究相对有限，主要集中在资源分布、增殖容量和资源评估方面。如徐炳庆等^[12]通过轮廓系数法对莱州湾三疣梭子蟹资源的时空分布特征进行了研究；栗小东等^[13]基于机器学习方法，对东海北部海域三疣梭子蟹的资源分布进行了预测；林群等^[14]构建了黄河口邻近海域的Ecopath模型，初步评估了三疣梭子蟹的增殖容量；俞存根^[4]对东海中北部重要的经济蟹类及其现存资源量和最大可持续产量(Maximum sustainable yield，MSY)进行了分析；Wang等^[15]基于传统Schaefer模型对三疣梭子蟹的最大可持续产量及资源开发情况开展了评估。以上在对三疣梭子蟹进行资源评估的研究中并未考虑模型中相关参数的不确定性。在资源评估过程中，如果不考虑不确定性，将对渔业资源评估结论、渔业生产和生产管理造成负面影响。因此，实现渔业科学管理和合理开发的前提是对资源量的正确评估^[16]。

本研究依据2001—2020年东海北部海域三疣梭子蟹的渔获量数据和单位捕捞努力量渔获量(Catch per unit effort, CPUE)数据，运用Schaefer剩余产量模型，将参数基准方案(参数的先验分布均为均匀分布)与敏感性分析方案(参数与基准方案相同，但采用不同先验分布)进行对比，对东海北部海域三疣梭子蟹的资源进行评估，并设定收获率作为三疣梭子蟹资源管理策略且分析管理策略可能存在的风险。研究结果有助于掌握当前三疣梭子蟹的渔业资源变动规律，并基于目前三疣梭子蟹资源状况对未来资源的动态变化趋势进行预估。

1 材料和方法 1.1 数据来源

东海北部海域三疣梭子蟹渔获量数据来源于《中国渔业统计年鉴》^[17]，时间段为2001—2020年(见图 1)。除了资源量以外，其他一些因素也会导致CPUE的变动，因此需要对CPUE进行标准化，以剔除这些因素的影响。本研究使用了近20年的渔业生产数据，其中包含了3种三疣梭子蟹生产船型，同时考虑主要环境因子变化的影响，因此在使用广义线性模型(Generalized linear model，GLM)对CPUE数据进行标准化时^[18]，考虑了年份、渔船类型和平均水温3个因素的影响。具体模型如下：

$\begin{array}{l} {\rm{lm(lnCPUE \sim factor(Year) + factor(Effort) + }}\\ {\rm{factor(Vessel) + factor(Temp), data = RData) }}。\end{array} $

(1)

式中: lm是R语言中广义线性模型的函数命令；ln CPUE代表对CPUE取对数；Year代表年份；Effort代表捕捞努力量；Vessel代表 3种不同的船型，分别是刺网、蟹笼、单拖网；Temp代表底层温度。

1.2 剩余产量模型

Schaefer非平衡剩余产量模型的离散时间形式如下:

$ B_t=B_{t-1}+r B_{t-1}\left(1-\frac{B_{t-1}}{K}\right)-C_{t-1} ; $

(2)

$ \begin{gathered} I_t=q B_t \mathrm{e}^{\varepsilon_t-\frac{\sigma^2}{2}}, \varepsilon_t \sim N\left(0, \sigma^2\right) 。\end{gathered} $

(3)

式中：B_t为t年的生物量；r为内禀自然增长率；K为环境容纳量；C_t-1为t-1年的渔获量；q为可捕系数。其中r、K、q为模型的参数。I_t为三疣梭子蟹资源丰度的相对指数，ε_t为误差项。

在进行资源评估的过程中，一般会假设剩余产量模型中初始年的生物量B₀为K，通过减少模型估算的参数数量来达到降低模型评估难度的目的^[19]。但这种假设成立的前提条件是能够确定渔业开始时前一年的资源量^[20]。由于中国对三疣梭子蟹的捕捞活动开始较早，因此将初始一年的资源量视为K并不合理。本研究的B₀值参考Wang等^[15]的研究结果。

1.3 似然函数

假设CPUE正比于生物量，且观测误差服从对数正态分布，则似然函数的表达式为:

$ \mathrm{L}(I \mid \theta)=\Pi_{t=2001}^{2020} \frac{1}{I_t \sigma \sqrt{2 \pi}} \exp \left(\frac{\left(\ln \left(I_t\right)-\ln \left(q B_t\right)\right)^2}{2 \sigma^2}\right) 。$

(4)

式中：L和ln均表示自然对数；I表示单位捕捞努力量渔获量；θ表示随机变量；Π表示连乘符号；q表示可捕系数；B_t表示t年的生物量；σ表示标准差。

由于本研究采用的渔获量和CPUE数据的时间序列为20年，这对标准差(σ)的估算来说存在一定的不稳定性，因此在本研究将其视为恒定值(σ=0.2)^[19]。

1.4 模型参数先验分布的设定

在贝叶斯理论中，参数的后验分布会受到先验分布设定的影响^[21]。因此，要对先验分布进行敏感性检验来验证样本数据提供的信息是否充足。具体方法是比较各参数在不同先验分布下得到的后验分布的差异大小，来判断样本数据提供的信息是否充足。若结果差异较大，则表明数据没有提供足够的信息^[18]。

因为较难搜集到三疣梭子蟹相关参数的先验信息，为了结果的合理性，故采用均匀分布作为模型参数先验概率分布的基准方案。此方案K的基准先验分布假设服从均匀分布(50×10⁴ t，550×10⁴ t)，下边界设置为50×10⁴ t，此值大于2001—2020年数据中的最大渔获量。为了尽量减少上边界值对后验概率的影响，上边界值被设定为550×10⁴ t^[22]。对于可捕系数q，获取其有信息的先验分布具有较大难度，在研究中一般将这种无信息的先验分布设置为均匀分布或对数均匀分布。为足以覆盖内禀增长率r的合理范围，r的上、下边界区间一般包括与评估种类相似的物种r的估计值^[23]。本研究为了验证样本数据是否提供了足够的信息，用敏感性分析方案来替代基准方案中的参数先验分布，比较两种情况下得到的参数后验分布差异大小。

根据三疣梭子蟹的相关研究^[15]，参数r、K、q基准方案和敏感性分析方案的先验分布设定见表 1。

1.5 模型后验分布的计算

贝叶斯方法不仅考虑了评估模型的不确定性，而且考虑了参数的不确定性，通过将经验信息与似然函数结合在一起，得到模型参数的后验分布，所以可以得到比传统分析方法更可靠的结果。因此，贝叶斯方法被认为是目前国际渔业资源研究领域较有代表性的方法，且适用于缺乏某些特定数据的情况下研究渔业资源动态^[24]。

本研究中模型参数r、K、q的后验概率分布是通过MCMC(Markov chain monte carlo algorithm)算法测算的。此算法共迭代50 000次，将前5 000次的结果舍弃，后每隔40次储存剩下的45 000次的运算结果。采用此算法计算时，各参数的初始值设定分别为：r=0.15，K=37 706.4，q=0.005^[15]。

1.6 生物学参考点估算

生物学参考点是从生物学角度设定的指标，包括目标参考点(Target reference points，TRP)和限制参考点(Limit reference point，LRP)。目标参考点是渔业管理为了持续获得某一渔获量所设定的指标，限制参考点指的是在渔业管理中为了避免出现种群的生物量过低导致种群不能被可持续利用的情况而设置的指标^[25]。本研究涉及到的生物参考点有F_MSY、B_MSY、F_0.1和MSY。F_MSY和B_MSY分别指渔业达到MSY水平时对应的捕捞死亡系数和生物量，F_0.1表示平衡渔获量和捕捞死亡系数关系曲线最大斜率的10%对应的捕捞死亡系数。通过计算可得：

$\begin{gathered} F_{\mathrm{MSY}}=\frac{r}{2}, \end{gathered} $

(5)

$ F_{0.1}=0.45 r, $

(6)

$ \mathrm{MSY}=\frac{r K}{4}, $

(7)

$ B_{\mathrm{MSY}}=\frac{K}{2} 。$

(8)

在本研究中，捕捞死亡系数的目标参考点F_tar设定为F_0.1，限制参考点F_lim设定为F_MSY；B_MSY设定为资源量的目标参考点B_tar，B_MSY/4设定为限制参考点B_lim。可以通过以上参考点来判断当前的渔业资源状况。当捕捞死亡系数大于F_lim时，表明该资源正在遭受过度捕捞(Overfishing)，相反则表明该资源未遭受过度捕捞；当资源量小于B_lim时，表明该资源已经遭到破坏，处于过度捕捞状态(Overfished)，相反则表明该资源未处于过度捕捞状态。

1.7 确定备选管理策略

本文采用不同的收获率(Harvest rate，即捕捞死亡率，规定了渔业年捕捞渔获量在资源量中的固定占比)来作为三疣梭子蟹资源管理策略。收获率分别设定为0.01、0.02、0.03、0.04、0.05、0.06、0.07和0.08。因此，未来t年的渔获量的计算公式如下：

$ C_t=h_i B_t \mathrm{e}^{\varepsilon}, \varepsilon \subset N\left(0, 1^2\right) 。$

(9)

式中：C_t为t年的渔获量; h_i为收获率; B_t为t年的生物量; ε为误差项。

1.8 管理策略实施效果模拟

资源评估过程的复杂性会导致评估结果存在较多不确定性，因此通过建立指标来对不同的管理策略进行风险评估是十分必要的。本研究假设实施管理策略的时间段为2021—2030年。建立的指标有：

(1) 管理结束时的生物量B₂₀₃₀，即2030年资源量的期望值。

(2) 管理策略实施后2030年渔获量的期望值C₂₀₃₀。

(3) 2030年三疣梭子蟹生物量与B_MSY之比的期望值(B₂₀₃₀/B_MSY)。

(4) 2030年三疣梭子蟹生物量的衰减率(Depletion)，即实施某一管理策略后，2030年三疣梭子蟹生物量与K的比例(B₂₀₃₀/K)。

(5) 2030年三疣梭子蟹生物量大于2020年三疣梭子蟹生物量的概率P(B₂₀₃₀＞B₂₀₂₀)，表示管理措施实施后资源较管理措施实施前资源恢复的概率。

(6) 2030年三疣梭子蟹生物量大于B_MSY的概率P(B₂₀₃₀＞B_MSY)，表示资源恢复到健康水平的概率。

(7) 2030年三疣梭子蟹生物量小于B_MSY/4的概率P(B₂₀₃₀ < B_MSY/4)，表示管理措施实施后资源崩溃的概率。

在1.5中，我们得出剩余产量模型参数r、K、q的后验概率分布后，可以通过反复的模拟计算得到指标概率分布，具体步骤见图 2。

2 结果 2.1 后验概率分布和期望值

图 3为不同方案模型的后验概率分布。基准方案假设r、K、q的先验概率均为均匀分布，但由图 3的曲线可知r、K、q的后验概率明显不同于假设的分布类型。同时，由两图对比可知这两种方案的参数后验分布非常相似，这表明模型参数后验分布的结果受本研究采用的渔业数据的影响远大于先验分布的影响，即数据提供了较丰富的信息。

模型估计的参数估计值见表 2，结果表明两方案下内禀增长率r的估计值分别为基准方案0.192，敏感性分析方案0.194，即后验分布均值无明显差异；基准方案和敏感性分析方案的环境最大容纳量K的估计值分别为503.38×10⁴ t和504.23×10⁴ t；两方案下可捕系数q的估计值相等。

2.2 三疣梭子蟹资源状况分析

基准方案和敏感性分析方案估算的三疣梭子蟹生物学和渔业学参考点见表 3。基准方案下估算的MSY为24.07×10⁴ t，MSY对应的生物量B_MSY为251.70×10⁴ t，此时的捕捞死亡系数为0.096。此方案评估的F_0.1为0.086。在敏感性分析方案下，估算的MSY为24.33×10⁴ t，B_MSY为252.12×10⁴ t，此时的捕捞死亡系数为0.097。此方案评估的F_0.1为0.087。由这两种方案估算得到的结果可知，F₂₀₂₀的均值小于F_MSY的均值且B₂₀₂₀的均值大于B_MSY的均值，这表明当前东海北部海域三疣梭子蟹资源未遭受过度捕捞且未处在过度捕捞的状态(见表 3)。

图 4的4个相位，红色区域代表已过度捕捞且正遭受过度捕捞，而绿色区域代表未出现过度捕捞，黄色区域代表已过度捕捞或正在过度捕捞^[26]。从图 4分析可知，2001—2020年，捕捞死亡系数基本低于F_MSY且生物量大都在B_MSY之上，这表明了近年来三疣梭子蟹资源状况良好，基本不存在过度捕捞的状况。

2.3 管理策略分析

基准方案和敏感性分析方案在不同的收获率下的管理决策和风险分析的指标结果见表 4。两种方案的各项指标值均比较接近，在基准方案下，随着收获率的增大，2030年渔获量期望值总体呈上升趋势，但当收获率从0.05增大到0.07时，2030年最大渔获量期望值的变化不显著。在敏感性分析方案下，收获率从0.01增大到0.06时，2030年最大渔获量期望值呈上升趋势，但随着收获率的继续增大，2030年最大渔获量期望值却逐渐减小。综合两种方案的结果可知，如果要获取2021—2030年最大持续产量，应将收获率控制在0.06左右，此时对应的资源量为23×10⁴ t左右。

当收获率为0.06时，风险分析表明资源状况并不十分乐观。在基准方案下，B₂₀₃₀/B_MSY的概率期望值为0.96，敏感性分析方案为0.99。对于B₂₀₃₀>B_MSY，基准方案和敏感性分析方案的概率期望值分别为0.53和0.56。对于B₂₀₃₀ <B_MSY/4，基准方案和敏感性分析方案的概率期望值分别为0.03和0.04。这表明该情况下，2030年的生物量与B_MSY非常接近且出现资源崩溃的可能，但生物量大于B_MSY的概率较小，此时的资源状况已出现衰退倾向，若继续保持0.06的收获率，将来资源崩溃的可能性会大大增加。

从渔业管理的角度出发，通过分析表 4可知，当收获率为0.04时，基准方案和敏感性分析方案中B₂₀₃₀>B₂₀₂₀的概率期望值均大于0.6，B₂₀₃₀>B_MSY的概率期望值均大于0.8，且B₂₀₃₀ <B_MSY/4的概率期望值较小，所以，将收获率控制在0.04左右是较为保守的决策，此时对应的持续渔获量在20×10⁴ t左右。此时资源量会维持在较高水平且不会出现资源崩溃的状况。

3 讨论与分析 3.1 先验和后验概率分布

贝叶斯方法是利用统计数据(如渔获量和CPUE等)与先验概率分布相结合，从而推断出后验概率分布的统计学方法^[27]。先验概率分布可以依据以往的研究结果，也可以参考研究者的主观判断。由于先验概率分布会对资源评估结果产生影响，选择了不合理的先验概率分布可能会导致错误的结果，所以对先验概率分布的推断目前仍有争议^[19]。因此，需要选择正确的先验概率分布，以避免后验概率分布结果完全被先验概率所主导^[28]。

在本研究中，三疣梭子蟹资源评估方面的参考资料较少，信息不足，只能依据已有的研究结果设定模型参数K、q、r的先验概率分布^[15]。因此，为了减少先验分布对后验分布的影响，在基准方案中，我们选择了范围较大的均匀分布，同时设置了敏感性分析方案进行比较分析。从后验分布来看，参数K、q、r与先验的均匀分布相比差异很大，说明本研究采用的数据提供了充足的信息；敏感性分析方案下的后验分布与基准方案下的后验分布十分相似，说明资源评估的结果受到模型参数先验概率分布的影响十分有限。

3.2 当前资源状况评价及未来管理策略分析

当前东海北部海域三疣梭子蟹资源状况相对乐观，在近年来并未遭受到过度捕捞且三疣梭子蟹的资源状况良好。本研究得出的捕捞死亡系数F较小，可能是由于三疣梭子蟹适应力强、存活率较高且增殖放流数量大，使得三疣梭子蟹资源量维持在较高水平，因此估算出的捕捞死亡系数F较小。本研究得出的三疣梭子蟹的MSY约为24×10⁴ t，这一结果不同于Wang等^[15]研究结果。Wang等^[15]研究显示，假设三疣梭子蟹的捕捞强度和放流数量保持在2015年的水平，三疣梭子蟹的MSY约为14.2×10⁴ t。造成结果不同的原因可能是采用了不同的评估方法或所用数据的时间跨度不同。Wang等^[15]使用的是传统的非平衡剩余产量模型，评估所使用的数据是2001—2015年间的渔业统计数据。当采用相同时间跨度的数据(2001—2015年)时，使用贝叶斯方法得到的MSY为19×10⁴ t，与Wang等^[15]14×10⁴ t的研究结果差异明显缩小。具体分析可得，2001—2015年的年产量均值约为10.4×10⁴ t，而采用2001—2020年的渔业统计数据算出的年产量均值约为12.1×10⁴ t。2015年的年产量约为18.5×10⁴ t，与本文评估得到的MSY值更接近。对于Wang等^[15]未使用的时间段即2016—2020年，东海北部海域三疣梭子蟹的年产量均值约为17.0×10⁴ t，明显大于2001—2015年的均值，且2001—2015年CPUE均值也偏小。对数据趋势分析可得，2013—2020年的CPUE均值明显大于2001—2012年，且2013年之后三疣梭子蟹年产量相较于之前明显增加，2015年之后的产量虽略有下降，但总体仍保持较高水平。因为本文采用的数据时间跨度为20年，故后期研究海域三疣梭子蟹统计数据的变化可能会对模型产生较大影响，导致评估结果(如MSY)数据变大。由此可见，采用不同时间段的数据做出的评估结果可能会产生差异，影响主要来自渔获量和捕捞努力量的年际变化趋势。因此，全面且准确的长时间序列数据对于评估MSY及相关的生物学参考点至关重要。

渔业资源管理的主要目标是通过合理的调控与科学的管理，从渔业资源中获得最佳的可持续利益。这就要求在对渔业管理策略进行选择时，应从多角度考虑，要在保证渔业资源可持续开发利用状态与获得高渔获量之间进行全面权衡。综合两种方案下的管理策略和风险指标分析结果可知，若要获取最大产量，应将收获率控制在0.06。但风险分析表明，此时的资源状况并不乐观。若将收获率设定为0.06开展捕捞活动，三疣梭子蟹资源会出现过度捕捞现象，甚至有可能会出现资源崩溃的情况。所以，将收获率设置为0.06对资源的可持续发展来说存在很大风险。

从资源保护的角度来看，将三疣梭子蟹资源的收获率控制在0.03时，资源量会维持在较高的水平且不会出现资源崩溃的状况。从经济收益的角度来看，将收获率控制在0.05时，基准方案和敏感性分析方案中2030年最大渔获量期望值与结果中最大值非常接近，且出现资源崩溃的可能性较低。但在大多数渔业已经出现衰退趋势的现状下，我们应采取使资源状况更安全的策略作为最佳管理策略。对利益及其存在的风险进行权衡后得出，0.04的收获率在未来10年的管理期内是较合理且保守的管理策略。

3.3 研究的不足和展望

本研究分基准和敏感性分析两种方案，使用了与贝叶斯方法相结合的Schaefer剩余产量模型，对东海北部海域三疣梭子蟹的资源状况进行了评估，并估算了三疣梭子蟹资源在2030年恢复到健康水平的概率和崩溃的概率。研究结果可为三疣梭子蟹资源养护和可持续利用提供有利支持，为制定科学合理的渔业管理策略提供参考。

但是，本研究仍然存在一些不足和需改善之处。第一，模型参数的不确定性，本研究假设了CPUE数据的标准差σ为0.2，该值的假设对评估结果造成的影响需要进一步研究与分析。第二，本研究采用的生产数据时间段为2001—2020年，因此模型参数的结果只能较客观地反映这一时期的渔业状态。由于三疣梭子蟹属于易受环境影响的短生命周期物种，且渔业统计年鉴中对该种的统计也存在误差，因此，基于2001—2020年间估算出的模型参数对其未来的渔业资源变化进行的预测会存在一定的不确定性。第三，三疣梭子蟹是东海北部海域重要的增殖放流种类之一，每年幼蟹放流数量达到几千万尾，而本研究并未将该影响纳入模型评估当中，这无疑会对资源评估结果产生一定的影响。今后的研究中，我们将继续积累更多的统计和调查数据，进一步完善评估模型，在评估过程中充分考虑环境因子和增殖放流等因素的影响，并尝试使用生态系统模型对三疣梭子蟹的MSY开展评估，为合理开发、有效保护三疣梭子蟹资源提供科学参考。