不同模态下气候变化对智利竹筴鱼补充量的影响

张畅，李纲，陈新军,2,3,4,5*

( 1.上海海洋大学 海洋科学学院，上海 201306；2.青岛海洋科学与技术试点国家实验室 海洋渔业科学与食物产出过程功能实验室，山东 青岛 266237；3.上海海洋大学 大洋渔业资源可持续开发教育部重点实验室，上海 201306；4.上海海洋大学 国家远洋渔业工程技术研究中心，上海 201306；5.上海海洋大学 农业农村部大洋渔业开发重点实验室，上海 201306)

1 引言

智利竹筴鱼（Trachurus murphyi）是东南太平洋重要的经济种类，于20世纪70年代开始大规模开发和利用，1995年达到历年最高产量，为496万t，其后多年渔业资源出现衰退，2013年产量仅35万t，为近30 a以来的最低产量，2019年捕捞产量又恢复到63.8万t[1]。作为沿海国家和远洋渔业国家的重要捕捞对象之一，其资源评估与管理工作受到世界各国和南太平洋区域渔业管理组织的极大关注，补充量作为研究种群资源量变动的一个重要因素，也是资源评估管理中不可忽略的重要因素。

通常，渔业资源的补充量除了取决于资源本身的再生产机制，还受到生态环境的制约[2]。在目前已有研究中，主要通过两个方面考虑环境对补充量的影响：一是使用一个或多个已知环境变量，构建短期的补充量预测模型[3]；二是考虑长期环境气候变化对补充量的影响，即考虑模态变动[4]。自然的模态变动可能是剧烈且突然的，与海洋气候联系密切，它有别与捕捞因素的影响，在渔业管理中有着十分重要的意义[5]。近年来，基于贝叶斯模型平均法构建的亲体-补充量模型在渔业中得到较好应用[6-7]，该方法不仅能分析各环境因子对资源补充量是否具有长期、稳健的解释能力，还能用于补充量预测。为此，本研究拟通过贝叶斯模型平均法，分析在不同模态中的环境因子对智利竹筴鱼的亲体补充量影响及其内在变化规律，为其资源评估和科学管理提供参考。

2 材料与方法

2.1 渔业与环境数据

（1）渔业数据。1971-2017年智利竹筴鱼亲体量与补充量源自南太平洋区域渔业管理组织（South Pacific Reginal Fisheries Management Organization, SPR FMO）评估报告[1,8]，以1龄鱼的资源量作为补充量，以参与繁殖的群体资源量作为亲体量。

（2）海洋环境数据。考虑智利竹筴鱼种群结构复杂，结合已有文献研究[9-10]，选取目前已知的产卵期（11 月至翌年2 月）最大产卵场（35°～45°S，90°～75°W）（图1）的海表面温度(Sea Surface Temperature, SST)、海表面盐度（Sea Surface Salinity, SSS）、海表面高度（Sea Surface Height, SSH）的平均值作为影响补充量的潜在影响因子，SST数据来源于https://oceanwatch.pifsc.noaa.gov/erddap/griddap/index.html?page=1&itemsPer-Page=1000，SSH与SSS数据来源于 http://apdrc.soest.hawaii.edu/las_ofes/v6/dataset?catitem=71。海洋尼诺指数 （Oceanic Niño Index, ONI） 以 Niño3.4 区 中 每 年11月至翌年2月的SSTA平均值来表示，时间范围为1971-2017年，数据来源于https://ggweather.com/ONI/oni.htm。太平洋年代际振荡(Pacific Decadal Oscillation,PDO）可直接造成太平洋及周边地区气候的年代际变化，对中上层鱼类资源造成影响，为此，本研究选取每年11月至翌年2月的平均太平洋年代际振荡，作为表征其资源量变化的环境因子，时间范围为1971-2017年，数据来源于http://research.jisao.washington.edu/pdo/PDO.latest。

图1 智利竹筴鱼产卵场Fig.1 The spawning ground of Chilean jack mackerel

2.2 模态检验

基于连续t检验的模态变动分析(Sequential T-test Analysis of Regime Shifts, STARS），进行补充量模态分析[11]。当模态改变时，可通过模态跃变指数（Regime Shift Index, RSI）反映，其计算公式为

式中，diff为两个连续且存在显著差异的格局模态平均值的差值；j为新模态的起始年份；L为所选取的截断长度；t为置信水平p下，2L-2自由度的t分布值；为方差；RSI为模态跃变指数；当模态向上转变，则=xi-， 当模态向下转变，则=-xi， 则=±diff；xi第i个变量；为L年内变量x的平均值；为新模态中变量x的平均值。根据相关文献和参数计算设定标准[11-12]，设置显著性检验置信水平p等于0.1，截断长度L等于20，huber权重参数H等于2。

2.3 补充量模型构建

选用Ricker模型构建基于环境因子的亲体-补充量模型[13]。其公式为

模型两边取对数后可转化为线性形式：

式中，R/S为繁殖率，又称为再生长率；α、β、β0为模型中的待定参数；R为补充量；S为亲体量；X表示影响因子；n为影响因子的数量。本次研究中变量因子包括亲体量、SST、SSS、SSH、ONI和 PDO，共 6个变量，可构建k=26个模型，所有模型构成模型集合M={M1，M2，···，Mk}。

2.4 贝叶斯模型平均法

2.4.1 模型参数估计

贝叶斯模型平均法（Bayesian Model Averaging，BMA）的基本原理为，根据已有的主观信息指定先验分布，解释变量的后验包含概率（Posterior Inclusion Probabilities，PIP）大小可作为选择解释变量的客观标准，以模型后验概率（Posterior Model Probabilities，PMP）为权重对2.3节中所有单项模型进行加权平均，形成包含所有信息的集合M，用于被解释变量的预测[14]。

为了不对模型的优劣加入主观色彩，模型先验设为均匀先验分布（即相同的模型先验概率），参数先验选择单位信息先验。利用马尔科夫链蒙特卡洛（Markov Chain Monte Carb，MCMC）方法在模型空间中抽样，计算迭代10 000次，抛弃前2 000次的计算结果。

模型后验概率PMP表达式为

各因子的后验包含概率PIP的表达式为

参数 β的后验分布为

式中，y为ln(R/S)；X表示影响因子；为模型Mγ的先验分布概率；为观测数据的概率总和。

2.4.2 补充量预测

贝叶斯模型平均法不仅可用于解释变量的选择，还可用于预测被解释变量。不同模态中，不同环境因子以不同的组合形式拟合多个亲体补充量模型，以PMP为权重对所有可能的模型进行加权平均，通过BMA综合多个模型进行集合预报，求得最后两年的预报结果[17]。为便于考虑模态变动对亲体-补充量的影响，将根据模态分析结果，设计多个基于BMA的亲体-补充量预测模型用于对比分析。

式中，yt为t年ln(R/S)预测值；yt,i为t年第i个模型的ln(R/S)预测值；为第i个模型的后验概率。

3 结果

3.1 模态检验

1971-2017 年，智利竹筴鱼补充量发生了显著的年际变化，大体上呈现先上升再降低，再上升降低的连续波动过程，具体表现为补充量在1971-1985年处于波动上升阶段，且在1985年达到最高值，在随后的几年里其补充量持续下降，1995年后再次上升，于2000年达到第2个高峰值，随后补充量下降，在2002-2015年维持在较低的范围内（图2）。STARS分析结果表明，1971-2017年间智利竹筴鱼资源补充量可分为5个模态：第1模态为1971-1980年，该模态中补充量处于上升阶段，补充量平均水平相对较低，多个国家和地区对该渔业进行了开发和利用；第2模态为1981-1990年，该模态中的补充量平均水平最高；第3模态为1991-2001年，该阶段补充量明显的起伏变化，但其平均值仍保持在较高的水平；第4模态为2002-2015年，该模态中的补充量平均水平最低；第5模态为2016-2017年，为新的模态阶段，与前一个模态相比，该阶段补充量水平有明显的上升，但该模态时间序列短，故在后续分析中不予讨论（图2，图3）。

图2 1971-2017年智利竹筴鱼补充量模态检测Fig.2 Regime shift detection in recruitment for different values of Chilean jack mackerel during 1971 to 2017

图3 1971-2017年智利竹筴鱼补充量跃变指数Fig.3 Regime shift index values of recruitment for Chilean jack mackerel during 1971 to 2017

3.2 气候变化特征

由图4可知，在第1模态中，厄尔尼诺累计值持续下降。在第2模态中，正负序列（不同气候事件）交替的频次高，变动趋势不稳定。在第3模态中，连续4 a为正数序列，厄尔尼诺累计值持续上升，经历连续两年减少，1997年迎来强厄尔尼诺，厄尔尼诺指数高达1.38，其后4 a为负数序列，累计值持续下降。在该模态中，厄尔尼诺指数累计值有明显的上升下降过程，变动趋势明显。在第4模态中，2002-2006年间整体趋势上升，2007-2014年整体趋势下降，但连续增长或下降的过程没有第3模态明显。分析发现，发生模态变化的年份都为厄尔尼诺年（图4）。

图4 1971-2015年间厄尔尼诺指数的年际变化Fig.4 Annual changes of El Niño index during 1971 to 2015

由图5可知，在第1模态中，PDO指数累计值持续下降。在第2模态中，1981-1988年PDO指数累计值持续上升，其后两年持续下降。在第3模态中，1991年PDO指数骤降至-1.78，整体趋势变动不明显，但增减变化的频率更高。在第4模态中，2003年PDO指数上升至1.86，其后小幅增长，2008-2014年整体趋势下降，2015年又上升至2.25。PDO指数的累计值在3个模态中都经历了增长再下降的过程，但在第2模态、第4模态中该过程更明显、更规律，在第3模态中增减变化的频率更高。分析发现，发生模态变动的年份前后，PDO指数都有明显的变化（图5）。

图5 1971-2015年间PDO指数的年际变化Fig.5 Annual changes of PDO index during 1971 to 2015

3.3 不同模态中影响因子的筛选分析

在不考虑模态变动的亲体-补充量关系中，亲体量的PIP为99%，在环境因子中，SSH的PIP最高，为76%。第1模态中，亲体量的PIP占34%，PIP最高的环境因子为SSS，其PIP值为46%；第2模态中，亲体量的 PIP为 86%，El Niño的 PIP值为 49%，为该模态中PIP值最高的环境因子；第3模态中，亲体量的PIP为100%，PIP最高的环境因子为PDO，其值为41%；第4模态中，亲体量的PIP为98%，PIP最高的环境因子为 El Niño，其值为 33%（表1）。

表1 不同影响因子在不同模态阶段的PIP（%）Table 1 PIP (%) of different impact factors on different regimes

3.4 预测模型

在不考虑模态变动的亲体-补充模型中，6个因子各自组合形成多个模型，表2列举PMP最高的前5个模型。模型1为基于SSH的亲体-补充量模型，后验概率最高，为36%。5个模型集合预报的累计后验概率达到了70%。

表2 基于BMA的PMP最高的前5个模型的后验概率（%）Table 2 Posterior distribution (%) of top five models with highest PMP of all models from BMA for no shift model

考虑模态变动对亲体-补充量的影响，依次去除第1模态、第2模态、第3模态，设计3个基于BMA的亲体补充量模型 (表3)：（1）方案 a，由第 2模态、第3模态、第4模态数据（1981-2015年）构建亲体补充量模型；（2）方案 b，基于第 3模态、第 4模态（1991-2015年）构建亲体补充量模型；（3）方案 c，由第 4模态（2002-2015年）构建亲体补充量模型。3个方案中，方案a中最佳模型的后验概率最高为24%，各个方案的累计后验概率基本相同(表3)。

表3 不同方案中基于BMA的PMP最高的前5个模型的后验概率（%）Table 3 Posterior distribution (%) of top five models with highest PMP of all models from BMA for different program

由图6可知，从预测结果看，BMA可通过概率分布定量评价预报的不确定性（图6）。不区分模态的预测结果最接近真实值，优于方案a、b、c，方案c最为保守，预测值明显小于其他方案。但其他方案的预测概率约为60%～80%，方案c的预测概率大于80%，优于其他方案。

4 分析与讨论

4.1 气候变化与补充量模态变动之间的联系

南太平洋智利竹筴鱼分布极其广泛，东起智利、秘鲁和厄瓜多尔沿岸，西至新西兰、澳大利亚沿海，在 30°～60°S、78°～177°W 广阔海域均有密集的智利竹筴鱼群体，该区域也称为“竹筴鱼带”[18]。该海域大尺度环流、中尺度海流以及各种海洋环境要素的变化活跃，智利竹筴鱼的补充量受到环境要素、气候事件的影响明显[19-20]。变动的生态环境也可能会引起竹筴鱼补充量的模态变动，从而给资源评估和管理带来不可忽略的影响[21]。

本研究基于STARS对智利竹筴鱼补充量进行模态分析，结果表明其存在多个模态，且在不同模态中，占主导地位的环境气候因子并不相同。在第1模态中，亲体量对补充量变动的解释概率明显偏低，各影响因子的解释概率也明显存在异常，可能的原因是该阶段渔业资源变化更多的受到捕捞因素的影响[22]，并不适用于分析补充量与环境的关系。其余模态与传统环境因子分析（不区分模态变动）相比，在环境因子的解释概率排序上有明显区别。在传统环境因子分析中，SSH占有极高的解释概率，而考虑模态变动的分析结果认为，各环境因子在不同的时期存在不同的排序。从统计结果分析，SSH的PIP明显偏大且优于其他环境因子，极少有环境因子的解释概率超过60%[23]，从目前已有的智利竹筴鱼与环境关系的研究结果分析，也更偏向于区分模态的分析结果[24-25]。传统环境因子分析往往将长时间序列的资源变动归因于一个关键因子上，而忽略了资源变动受到多个环境因子的共同影响，在不同时期占主导地位的环境因子也会发生改变。要深入了解渔业资源对环境气候的响应以及鱼类动态变动及生态系统调节机制等内容，模态分析是不可忽略的。

在不同环境时期，智利竹筴鱼的补充量水平会发生改变，基于这些不同时期的平均补充量水平来拟合亲体-补充量关系和制定管理方案，可能不利于渔业的可持续发展[26]。在智利竹筴鱼的评估报告[1]中，选取两个基于不同时间段亲体-补充量关系（1970-2015年、2000-2015年），并分别放入渔业资源评估模型，其结果存在较大的区别。研究认为，基于2000-2015年数据的分析结果更为保守，可用于最终的资源评估与管理。此外，SPRFMO评估报告[1]中选取的时间段（2000-2015年）与本研究中的第4模态相近，基于第4模态的预测结果也明显比其余方案更保守（图6）。但从预测结果看，4种预测方案中，不考虑模态变动的方案最接近真实值，基于第4模态的结果偏差最大，其原因可能是用于预测的2016年和2017年补充量水平高于第4模态补充量水平，与2015年以前有较大的差异（图2），因此只基于第4模态（低水平）进行预测其结果自然会偏小。目前，在SPRFMO智利竹筴鱼评估报告中，数据段的选取更多的依赖于人为经验，但本文基于模态检验的科学方法将智利竹筴鱼划分为多个模态进行分析，能为数据段的选取及资源评估提供了科学依据。

图6 基于BMA的ln(R/S)预测概率Fig.6 Predictive probabilities for ln(R/S) based on BMA

4.2 厄尔尼诺和PDO事件对智利竹筴鱼补充量的影响

在第2模态、第4模态中，厄尔尼诺皆为解释概率最高的环境因子，其指数的正负序列（不同气候事件）交替频次高，持续处于某一气候事件的时间短，与第3模态中的表现有明显的区别。这可能是由于不同气候事件高频的交替发生，这种不稳定的气候环境对智利竹筴鱼资源补充量产生了更大的影响。研究认为，智利竹筴鱼在面对不同气候事件时，存在不同的应对机制。以厄尔尼诺和拉尼娜为例，在厄尔尼诺期间，海表面温度、海表面高度、风速、盐度等众多环境要素都会发生改变，智利竹筴鱼仔稚鱼分布、洄游路径等多方面产生变化[27]。同时，由于捕食者减少、不同性质海水扩展等原因，智利竹筴鱼栖息范围会有所扩展[9]。Arcos等[28]研究认为，厄尔尼诺对智利竹筴鱼仔稚鱼肥育场的影响往往会持续3～4 a的时间。在拉尼娜期间，气候条件与厄尔尼诺存在显著不同，海水的降温可能导致智利附近海域的幼鱼进入厄瓜多尔和智利的专属经济区[29]。

PDO是一种反映从年际到年代际的气候变率强信号，是太平洋重要的气候系统之一，对太平洋乃至全球的大气环流有着重要影响，也是气候预测的关键因子[30-31]。除了对各地区气候的影响，PDO对海洋鱼类种群的兴衰也起着调节作用[32]。PDO与许多中上层鱼类的补充量有着直接的联系，PDO正负相位时期、冷暖时期对应的补充量都有明显的变化[33]。Dioses等[34]研究认为，年际间PDO变化与智利竹筴鱼的纬度隔离有一定关系。分析认为，1977-1999年为PDO暖期，1999-2015年为PDO冷期[35]。冷暖期交替发生在本研究中的第3模态时期，这也解释了该模态中PDO年际间变动更为不稳定，而不稳定的年际间PDO变化对智利竹筴鱼补充量产生了更大的影响。

本研究认为，在发生模态变动年份前后，均发生了厄尔尼诺事件且PDO指数有明显的跃变，这表明智利竹筴鱼补充量与厄尔尼诺和PDO气候事件存在紧密联系，PDO年代际冷暖期交替与厄尔尼诺现象等可能是诱发智利竹筴鱼补充量发生模态转变的重要因素。Espino[36]研究也表明，PDO、厄尔尼诺、智利竹筴鱼三者间的变动存在紧密联系。一方面，不同尺度间的气候变化存在交互作用，年代际PDO变化会影响厄尔尼诺现象，而物理环境中这些低频的变化可能会导致海洋生态系统的改变[37-38]，从而导致智利竹筴鱼补充量的模态改变。另一方面，PDO和厄尔尼诺会以不同的方式和不同的量级对渔业资源造成影响[14]。在今后的研究中，需要结合物理海洋、气候学等多学科，加强了解智利竹筴鱼补充量对气候环境变化的响应机制和机理。

4.3 其他环境因素对智利竹筴鱼补充量的影响

海表面高度是用于检测和表征中尺度涡流的重要物理量[39-40]。涡流系统对智利海域叶绿素的聚集起着重要作用，冬季高浓度叶绿素海域多数是因此形成[41]。同时，涡流和上升流之间的耦合使得高初级生产力区域扩张，持续时间延长，浮游生物随之聚集，为智利竹筴鱼仔鱼提供连续的食物来源，有利于其生长[42-43]。此外，涡流的垂直影响深度完全覆盖仔鱼生活水层，部分仔鱼被夹带入涡流内会被运到外海，在适当条件下可能会促进产生新的仔、稚鱼栖息地[44-45]。在早期生活史阶段，涡流等中尺度海流变化过程对智利竹筴鱼空间分布、生长以及资源补充量等有着不可忽略的影响。

在本研究中，温盐的PIP基本维持在20%～30%，这表明两者对智利竹筴鱼补充量也具有一定的影响。东南太平洋海域拥有众多不同温盐特质的水团，占主导地位的水团改变及强弱改变都会导致竹筴鱼分布情况的改变[46]。当高温（温度大于18.5℃）、高盐（盐度大于34.9）亚热带表层水团南移并覆盖近岸海域时，鱼群分布更为集中，渔获量增加；当低温（约12℃）、低盐（约34.25）的亚南极中层水团强于赤道表层水时，智利竹筴鱼分布在更深的水层中，渔获量明显减少。秘鲁海域的智利竹筴鱼常聚集在表温为11～28℃的水域，智利和公海海域智利竹筴鱼聚集在表温为9～20℃的水域, 当海表温上升时，智利竹筴鱼适宜栖息地有明显的南移趋势。与其他鱼类不同，盐度对智利竹筴鱼的影响更为直接[47]，张敏等[29]研究认为，较稳定的盐度将有利于智利竹筴鱼集群的稳定。

模态检验方法的参数设置不仅需要遵循模型使用的前提条件，同时还需考虑智利竹筴鱼渔业资源本身变动规律及其特性，任意降低截断时间长度可能会出现更多的模态[48]，时间序列过短的模态不利于渔业资源评估和管理，选择一个恰当合理的时间序列长度在不同模态的渔业资源评估是极为重要的。气候环境变化对智利竹筴鱼早期生活史阶段的影响是不可忽略的，合理的划分模态，有助于探究不同阶段的环境与补充量的影响模式。其次，利用模态变动的方法将长期气候变化对补充量的影响纳入考虑，有助于提高对鱼类动态及其生态系统调节机制的理解，制定更合理的管理策略[49]。同时，贝叶斯模型平均法将补充量影响机制归因于多个环境因子共同作用的结果，分析每个因子不同的解释能力，还能以概率分布的形式描述补充量预报的准确性。模态分析与贝叶斯模型平均法的结合适用于多信息集合，可在今后的研究中对更多的环境气候因素进行分析，促进物理海洋、气候、生物等多学科交叉，建立起科学的智利竹筴鱼补充量预报模型，为后续的渔业资源评估和管理提供基础。