气候变化对人类社会和自然生态系统的威胁是目前最重要的国际问题之一，气候变化通过对地球物理系统的自然影响来作用于其中的各种生物。作为对正在发生的气候变化的响应，许多陆地、淡水和海洋物种已改变了其分布范围、季节性活动、迁徙模式、丰度以及物种间的相互作用^[1]。研究表明大西洋和太平洋的主要商业鱼类捕获量变动都与长期的气候波动有关^[2]。

智利竹筴鱼(Trachurus murphyi)是一种分布于南太平洋亚热带水域的上层鱼类，从南美洲西海岸的秘鲁、智利一直到新西兰、澳大利亚外海均有广泛分布，是秘鲁、智利两国沿岸围网渔业，也是我国远洋渔业大型拖网加工船的主要目标鱼种之一^[3]。智利竹筴鱼渔获量在20世纪70年代开始稳步增加，到1995年达到高峰后便迅速下降，之后一直在低位徘徊，其中有全球气候变化在内的多种原因^[2]。目前已有一些学者对智利竹筴鱼资源量受气候变化或海洋环境因素的影响进行了研究，如Espíndola等^[4]在运用Ricker模型对智利竹筴鱼进行资源评估时发现，智利竹筴鱼补充量受海表面温度的影响; Li等^[5]以及牛明香等^[6]进行智利竹筴鱼渔获量标准化时得出智利竹筴鱼资源丰度受海表面温度、拉尼娜和厄尔尼诺事件以及叶绿素浓度等的影响; Arcos等^[7]的研究也说明了智利竹筴鱼资源分布受厄尔尼诺事件的影响等。已有研究中智利竹筴鱼渔业受气候变化影响涉及的因素以海洋环境各要素为主，对区域以及全球气候变量的讨论较少，原因在于气候变化往往持续较长一段时期，通常为数十年或更长时间^[8], 而目前的研究多以10年以内的短期研究为主，短期研究在揭示气候变化对渔业的影响上存在一定的局限性。本研究采用低频气候变化参数长时间序列数据，对其进行相关性分析，并在此基础上利用神经网络模型构建东南太平洋智利竹筴鱼渔获量预测模型，用以探讨气候变化对智利竹筴鱼渔获量的长期影响，为东南太平洋智利竹筴鱼的可持续开发提供科学的参考依据。

采用低频气候变化参数: 厄尔尼诺相关指标(Niño1+2、Niño3、Niño4以及Niño3.4), 南方涛动指数(Southern Oscillation Index, SOI), 北大西洋涛动(North Atlantic Oscillation, NAO), 太平洋年代际涛动(Pacific Decadal Oscillation, PDO), 北太平洋指数(North Pacific Index, NPI)以及全球海气温度异常指标(dT)作为气候变化表征因子。其中，dT为英国气象局哈德利中心的观测数据集Hadcrut4数据(https://www.metoffice.gov.uk/hadobs/ hadcrut4/), 该数据是由分布在全球50多个国家的观测站观测的、包含陆地表面气温和海表面温度的100个数据集融合而来的网格数据; 其余数据来自美国国家海洋和大气管理局(https://www. esrl.noaa.gov/psd/gcos_wgsp/Timeseries/, 各指标含义可以参考上述网站), 各指标数据均为1900— 2016年月度数据。

东南太平洋智利竹筴鱼总渔获量(catch)为南太平洋区域渔业管理组织(South Pacific Regional Fisheries Management Organization, SPRFMO)管辖范围内的沿岸围网渔业以及公海拖网渔业所获全部渔获量(http://www.sprfmo.int/assets/SC5-2017/ Annex-9-2017-Assessment-input-and-results.pdf), 为1970—2016年年度数据。

变量间的相关性分析采用Spearman秩相关系数来衡量，这是一种无参数检验方法，不仅可度量变量间的简单线性关系，也可度量复杂的单调关系的强弱，因为渔获量与气候变化表征因子之间的关系往往非简单的线性关系，故本研究选用Spearman秩相关系数来度量变量间的相关关系，公式^[9]为:

式中，${d_i}$为两个按大小次序标记好等级的序列中每个序列对等级编号之差，$n$为样本容量。

一般$|{r_s}| \le {\rm{1}}$, $|{r_s}|$越接近1表明两变量间相关性越强，其显著性可由t统计量进行判断^[9]。

互相关分析可以描述随机变量X(·)与Y(·)在任意两个不同时刻的取值之间的相关程度^[10], 反映了两个变量在不同的相对位置上的相似程度，其表达式为:

式中，t为延迟量，X(·)与Y(·)为时间序列，E(·)为期望值函数。R(t)描述了在延迟变量t处、时间序列X(s+t)与Y(s)的相似程度，如果二者相关程度大，则R(t)较大。

鱼群对气候变化的响应往往具有一定的滞后性^[11-12], 采用互相关分析可以帮助判断气候变化表征因子对于竹筴鱼渔获量影响最大时的滞后年限。

采用BP神经网络模型对东南太平洋智利竹筴鱼渔获量进行预测。对于一个输入层、隐含层以及输出层分别含有n、m和l个分量的三层BP神经网络，其模型^[13]如下:

${y_j} = f\left( {\sum\limits_{i = 1}^n {{v_{ij}}{x_i}} } \right)\;\;\;(j = 1,2, \cdots ,m)$

${o_k} = g\left( {\sum\limits_{j = 1}^m {{w_{jk}}{y_j}} } \right)\;\;(k = 1,2, \cdots ,l)$

$e = \frac{1}{2}\sum\limits_{k = 1}^l {{{({d_k} - {o_k})}^2}} \;$

式中，x_j、y_j、O_k 分别为输入层输入分量、隐含层输出分量以及输出层输出分量，v_ij、w_jk分别为输入层分量对隐含层分量以及隐含层分量对输出层分量权重值，e为网络输出${o_k}$与期望输出${d_k}$的误差。f(·)与g(·)为转换函数。本研究设定隐含层转换函数为对数Sigmoid函数$f(x) = 1/(1 + {e^{ - x}})$, 输出层转换函数为线性函数$g(x) = x$。

隐含层节点数可依照经验公式^[13]来设置:

式中，m为隐含层节点数，n为输入层节点数，l为输出层节点数，α为1~10之间的常数。

由模型可知，误差最终是各层权重的函数，神经网络的训练即是通过误差方向传导机制不停地调整各层权重值(${w_{jk}},\;{v_{ij}}$), 以使误差e最小化。

对于神经网络模型进行性能评估可以采用效率系数^[14] (coefficient of efficiency, E), 效率系数反映模型预测值与实际观测值的平均值之间的接近程度，一般E≤1, 其值越接近1表示拟合值与真实值的一致性程度越高，模型越优良。

样本量为n, 输出样本为${o_i}$, 输出样本均值为$\bar o$, 模型计算的拟合值为${\hat o_i}$, 则效率系数为:

变量的敏感性分析是一种定量描述模型输入变量对输出变量重要性程度的方法^[15], 本研究采用Garson算法的改进算法^[15-16], 从神经元间的连接权值入手进行变量敏感性分析，若输出层只有一个输出变量，则第$i\;\;(i = 1,2, \cdots ,n)$个输入变量对此输出变量的敏感系数值为:

式中，${v_{ij}}$是输入层第i个变量对隐含层第j个神经元的权重，${w_j}$是隐含层第j$\;(i = 1,2, \cdots ,n;j = 1,$$2, \cdots ,m)$个神经元对输出层变量的权重。

由公式(5)可以看出，各个输入变量对输出变量重要性程度衡量的敏感系数值之和为1, 即$\sum\limits_{i = 1}^n {{Q_i} = 1} $。从而可以依据敏感系数值的大小来区别输入变量对指定输出变量影响程度的大小。

气候变化表征因子为${x_{ij}}$, 表示指标在第i年第j月的指标值，运用公式${x_i} = \frac{1}{{12}}\sum\limits_{j = 1}^{12} {{x_{ij}}} $得到该指标第i年的年平均数据。9个气候变化表征因子的时间序列图如图1所示。

可以看出，在1个多世纪的时间内，dT数据呈现出一种震荡向上的趋势; PDO数据呈现出一种长寿命周期(40年左右)的冷暖相位交替变化，而其余各个因子数据变化周期稍短。

东南太平洋智利竹筴鱼渔获量数据显示(图2), 该鱼种渔获量自1970年开始直线上升，直至1995年达到历史最高产量495.5万t, 随后其渔获量断崖式下降，直到2000年左右有所回升，但2008年后渔获量又呈下降态势。

对9个气候变化表征因子以及智利竹筴鱼渔获量数据运用线性变换进行标准化，公式为:

式中，${\bar x_i}$为变换后的因子数据，${x_{\max }}$为因子数据最大值，${x_{\min }}$为因子数据最小值，${x_{{\rm{mid}}}} = \frac{{{x_{\max }} + {x_{\min }}}}{2}$为因子数据中间值，变换后各指标值均统一在 [–1, 1], 标准化后的因子仍沿用原记号。

图1

图1 　气候变化表征因子的时间序列图 Fig. 1 　Time series of characterization factors of climate change

图2

图2 　东南太平洋智利竹筴鱼渔获量时间序列图 Fig. 2 　Time series of Trachurus murphyi catch in the Southeast Pacific

首先对气候变化表征因子Niño1+2、Niño4、Niño3.4、Niño3、SOI、NAO、PDO、NPI以及dT的样本数据进行相关性分析，以剔除其中相关性较大的因子，形成相对独立的因子指标体系。采用公式(1), 计算各个气候变化表征因子1900—2016年间的相关系数(表1)。

表1

表1 　气候变化表征因子Spearman秩相关系数表 Tab. 1 Spearman rank correlation coefficient of characterization factors of climate change 表征因子characterization factor Niño1+2 Niño4 Niño3.4 Niño3 NAO SOI dT PDO NPI Niño1+2 1.000 Niño4 0.467 1.000 Niño3.4 0.711 0.934 1.000 Niño3 0.892 0.783 0.943 1.000 NAO 0.218** 0.182** 0.241** 0.256** 1.000 SOI −0.583 −0.876 −0.909 −0.812 −0.289** 1.000 dT 0.153** 0.411 0.285* 0.260* −0.087* −0.143 1.000 PDO 0.435 0.589 0.580 0.528 0.154** −0.626 0.173** 1.000 NPI −0.245 −0.286 −0.257 −0.238** 0.048 0.335 −0.091 −0.655 1.000 注：无上标表示在5%显著性水平下通过检验; *表示在10%显著性水平下通过检验; **表示在10%显著性水平上未通过检验. Note: No superscript means passing the test at the 5% significance level; * means passing the test at the 10% significance level; ** means that it fails the test at the 10% significance level.

表1 　气候变化表征因子Spearman秩相关系数表 Tab. 1 Spearman rank correlation coefficient of characterization factors of climate change

从9个因子的Spearman相关系数来看，因子Niño1+2、Niño4、Niño3.4以及Niño3有相对较强的相关关系，且SOI与Niño4、Niño3.4以及Niño3之间有高度相关关系。从因子独立性的角度来看，可以从9个因子中除去因子Niño4、Niño3.4以及Niño3, 得到相对独立的气候变化表征因子集{Niño1+2、NAO、SOI、dT、PDO、NPI}。

采用公式(2)对智利竹筴鱼渔获量catch与因子集{Niño1+2、NAO、SOI、dT、PDO、NPI}中因子进行互相关分析，利用各因子1970—2016年年平均值数据，考虑气候变化因子的15年滞后(t=−1~−15), 比较不同滞后阶数时各因子与catch的互相关系数绝对值，互相关系数绝对值最大时对应的滞后阶数为气候变化因子的最佳滞后阶数(表2)。从表2可以看出，各个因子与智利竹筴鱼catch的相关程度各不相同，相关程度最大的是取3年滞后的因子PDO, 最小的是9年滞后的因子Niño1+2。

表2

表2 　各气候变化表征因子的最佳滞后阶数及其与渔获量的互相关系数值 Tab. 2 　The best lag order of climate change characterization factors and cross-correlation coefficient between climate change characterization factors and catch 表征因子characterization  factor 最佳滞后阶数the best lag  order 互相关系数值cross-correlation  coefficient PDO 3 0.501 dT 12 −0.468 SOI 2 −0.461 NAO 4 0.371 NPI 3 −0.281 Niño1+2 9 0.134

表2 　各气候变化表征因子的最佳滞后阶数及其与渔获量的互相关系数值 Tab. 2 　The best lag order of climate change characterization factors and cross-correlation coefficient between climate change characterization factors and catch

将2.3得到的6个气候变化表征因子滞后数据作为输入变量，智利竹筴鱼渔获量数据作为输出变量，建立一个三层BP神经网络模型。依照公式(3)计算得到隐含层节点数m为3~13, 因此制定了11套模型运行方案($m = 3,\;4,\; \cdots ,\;13$), 每套方案训练30次，共运行330次。

采用MATLAB (R2013b)软件的神经网络拟合工具建立BP神经网络，并实现其功能: 将47组(1970—2016年)样本按照分层抽样的方法分为训练样本、验证样本以及测试样本3部分，其中2008—2016年样本为测试样本; 1970—2007年样本每4年为一组，每组第一个样本设定为验证样本，其余3个样本为训练样本; 训练样本、验证样本以及测试样本三部分比例为28 : 10 : 9。选择LM (levenberg-Marquardt)算法进行样本训练，LM算法是一种二阶非线性优化算法，它具有非常快速的收敛性。设置网络最大训练次数为1000次，但为了避免过度拟合，并使训练更有效率，在模型运行过程中记录每个训练周期验证样本数据集的识别率，当其不再上升的时候，就停止训练。

为了使挑选出来的神经网络模型有较强的泛化能力，更具通用性，文中使用验证数据集效率系数E₁和测试数据集效率系数E₂[公式(4)]来选择最优模型。首先，在每一套运行方案的30次运行中依照E₁+E₂最大原则挑选出每套运行方案的最佳模型; 然后，比较11套方案的最佳运行模型的E₁+E₂值(表3), 其中最大值对应的方案为智利竹筴鱼渔获量预测最优BP神经网络模型。

可以看出，隐含层神经元m=6时的最佳运行模型为智利竹筴鱼渔获量预测最优BP神经网络模型。

将6个气候变化表征因子以及智利竹筴鱼渔获量数据值输入智利竹筴鱼渔获量预测最优BP神经网络模型，得到渔获量拟合值。比较模型拟合值与观测值(图3, 其中2008—2016年数据为测试样本), 可以看出，渔获量拟合值与观测值变化趋势基本一致，两个序列的线性相关系数为0.745; 对渔获量观测值与拟合值进行线性回归，用x表示渔获量拟合值，y表示渔获量观测值，得到回归方程为: y=−785.9+ 1.459x (图4直线), 其中方程的自变量显著性检验t值为15.2 (其伴随概率P值为0), 方程拟合优度R²=0.554, 这表明回归方程显著，最优BP神经网络模型拟合效果良好。

表3

表3 　BP神经网络模型的效率系数 Tab. 3 　Efficiency coefficient of BP neural network model 方案 scheme (m=) E 1 E 2 E 1+E2 3 0.5109 0.6077 1.1185 4 0.4241 0.7524 1.1766 5 0.4925 0.6510 1.1434 6 0.4281 0.8280 1.2561 7 0.2401 0.6766 0.9167 8 0.3189 0.6628 0.9817 9 0.3291 0.6272 0.9563 10 0.3649 0.5756 0.9406 11 0.0924 0.7042 0.7966 12 0.1959 0.8793 1.0752 13 0.2994 0.4713 0.7707

表3 　BP神经网络模型的效率系数 Tab. 3 　Efficiency coefficient of BP neural network model

图3

图3 　智利竹筴鱼渔获量观测值与拟合值对比图 Fig. 3 　Comparison between the observed catch and the fitted value of Trachurus murphyi

图4

图4 　智利竹筴鱼渔获量观测值与拟合值线性关系图 Fig. 4 　Linear relationship between the observed catches and fitted values of Trachurus murphyi

记录最优BP神经网络模型(m=6)中输入层到隐含层之间的权值矩阵为$V = {({v_{ij}})_{6 \times {\rm{6}}}}$以及隐含层到输出层之间的权值向量$W = {({w_j})_{{\rm{1}} \times {\rm{6}}}}$, 其中${v_{ij}}$为第i个输入向量对第j个中间变量的权重，${w_j}$表示隐含层第j$(i,j = 1,\;2,\; \cdots ,\;6)$个神经元对输出向量的权向量，最优模型的各级权重见表4。

按照公式(5)计算得到模型6个输入变量Niño1+2、NAO、SOI、dT、PDO、NPI对输出变量Catch的敏感性系数分别为0.14、0.23、0.08、0.11、0.23、0.21。从而，最优模型6个气候变化表征因子对智利竹筴鱼渔获量的影响大小从高到低依次为NAO、PDO、NPI、Niño1+2、dT以及SOI, 可以看出前3个因子的敏感性系数值非常接近，而后2个与前3个有一定的差距。

表4

表4 　气候变化表征因子对智利竹筴鱼渔获量影响的最优模型权重值 Tab. 4 　Weight values of optimal model of the impact of climate change characterization factors on Trachurus murphyi catch 表征因子characterization factor i j 1 2 3 4 5 6 v ij 1 0.53 −0.66 −1.25 −0.07 0.53 0.74 2 −1.28 0.19 0.61 1.09 1.18 0.36 3 0.593 −0.51 −0.17 0.03 0.43 0.60 4 −0.83 1.12 −0.36 −0.14 −0.52 0.54 5 0.60 0.20 0.60 0.92 −0.67 −0.72 6 0.32 0.89 0.82 0.44 −0.85 −1.09 w j −0.08 −0.09 0.07 −0.29 0.71 0.10

表4 　气候变化表征因子对智利竹筴鱼渔获量影响的最优模型权重值 Tab. 4 　Weight values of optimal model of the impact of climate change characterization factors on Trachurus murphyi catch

在渔业问题的研究中，气候变化表征因子越来越受重视，有学者甚至提出利用气候变化表征因子帮助管理渔业、对行业进行长期规划^[17]。对于本研究提取出的与智利竹筴鱼相关的3个主要相关因子NAO、PDO以及NPI, 已有众多研究表明它们对各种海洋生物的影响，如Beaugrand^[18]的研究发现大西洋内多种桡足类动物的丰度及多样性与NAO相关; Southward等^[19]的研究发现早在16世纪，英国德文郡和康沃尔郡的鲱(Clupea harengus L.)和沙丁鱼(Sardina pilchardus)渔业就受NAO影响，而且直到19世纪鲱和沙丁鱼丰度依然受其影响; Cabrero等^[20]研究发现在伊比利亚大西洋大陆架上NAO是影响沙丁鱼(Sardina pilchardus)补充量的主要因素; Sanchez-Rubio等^[21]的研究发现墨西哥湾北部的大鳞油鲱(Brevoortia patronus)上岸量与NAO呈负相关，Gamito等^[22]的研究发现葡萄牙中部塞图巴尔渔港欧洲沙丁鱼(Sardina pilchardus)的单位捕捞努力量渔获量(landings per unit of effort, LPUE)与NAO呈负相关关系; 而Sanchez-Rubio等^[23]的另一项研究发现NAO以及降水量等因素解释了墨西哥湾中北部蓝蟹(Callinectes sapidus) 约25%丰度的变化。

关于PDO对渔业的影响，有研究表明，2000年以来在厄瓜多尔水域捕获的鳀(Engraulis ringens)渔获量与PDO有较高的相关性^[17]; 此外，研究也表明在太平洋西岸，阿拉斯加与美国西海岸太平洋鲑在20世纪后70年渔获量发生的更替性变化与PDO正负相位转换有关^[24]; 与此同时，有研究指出在太平洋东岸，PDO是影响日本沙丁鱼(Sardinops melanostictus)和鲐(Scomber japonicas)补充量的一个重要因素^[25]; Li等^[26]在探讨气候对黄海鲱(Clupea pallasii)种群影响时发现黄海鲱的种群丰度与PDO等指标密切相关;  Castillo等^[27]的研究发现秘鲁外海鳀(Engraulis ringens)的分布可能与PDO等气候变化有关; Lennert-Cody等^[28]研究发现赤道东太平洋西部的小型和中型丝鲨(Carcharhinus falciformis)丰度变化与PDO相关; Shanks^[29]的研究发现2007年俄勒冈州库斯湾(Coos Bay)的太平洋大蟹(Cancer magister)大眼幼体渔获量的阶跃式变化与PDO相关; Zwolinski等^[30]研究发现加利福尼亚海流中的太平洋沙丁鱼(Sardinops sagax)的种群数量变化与PDO相关，并依此提出了一个基于季节性PDO指数的沙丁鱼补充量预测模型。关于NPI对渔业的影响，Goyert等^[31]利用信息论的方法评估了NPI等气候变化指标对阿拉斯加海鸟数量变化的影响，指出NPI与海鸟种群数量变化呈负相关关系; 也有研究发现，亚北极分水岭红鲑(Oncorhynchus nerka)的迁移持续时间与NPI的变化有较强相关性^[32]。可见，NAO、PDO、NPI对多种海洋生物包括浮游生物、鱼类以及鸟类的种群分布、季节性活动、丰度等都产生了不可忽视的影响。

另外，NAO、PDO等气候变化因子作为北半球气候变化的显著模态不仅影响着北半球的气候^[33], 还对南半球有着一定的影响，如Wu等^[34]利用观测资料和海气耦合模式研究了北太平洋的全球气候影响，发现北太平洋西部的海气耦合相互作用可能通过NAO等大气遥相关来调节热带印度洋、副热带西太平洋以及南太平洋的气候，即北太平洋可能是全球十年气候变化的重要调节器; 也有研究发现，在相隔数千公里的大型海洋生态系统中，主要鱼类种群在年代际尺度上经常出现同步波动，如黑潮和洪堡海流生态系统之间的物理和生物过程的同步性^[35], 科学家推测这种同步性是由南北太平洋海盆尺度耦合海气环流的变化引起的，并指出NAO似乎是这种物理生物同步性的同步剂。到目前为止，气候变化表征因子之间的相互作用以及气候变化对海洋生物影响的机理尚未完全清晰。本研究中，作为北半球气候变化显著模态的NAO、PDO以及NPI等对东南太平洋智利竹筴鱼渔获量的影响机理也有待进一步探讨。

本研究探讨东南太平洋智利竹筴鱼渔获量受气候变化的长期影响，将气候变化表征因子作为唯一的影响因素进行研究，这有别于目前已有的诸多短期研究，短期研究中往往将气候变化变量当作影响因素之一对智利竹筴鱼资源量进行探讨，如Espíndola等^[4]在运用Ricker模型对智利竹筴鱼进行资源评估时引入了海表面温度，牛明香等^[6]进行智利竹筴鱼渔获量标准化时在广义可加模型(GAM)中引入海表温度、叶绿素浓度和海表温度梯度，Li等^[5]在进行智利竹筴鱼渔获量标准化时引入了拉尼娜和厄尔尼诺事件等。仅以气候变化表征因子作为影响因素来进行研究，凸显了智利竹筴鱼渔获量变化的气候原因。另外，短期研究受数据时间序列长度的限制，不适合选择变化周期较长的低频气候变化参数作为影响因素，而常选择海洋环境变量如海表面温度等进行研究，或者选择某个低频气候变化参数出现模式转换的时间区间进行研究，如Arcos等^[7]研究了1997—1998强厄尔尼诺事件的发生对智利竹筴鱼资源分布的影响，这有别于采用低频气候变化参数长时间序列作为影响因子的渔业研究。采用低频气候变化参数长时间序列对东南太平洋智利竹筴鱼渔获量进行研究，可以从全球气候变化的角度来分析智利竹筴鱼渔业资源受气候变化的长期影响，为智利竹筴鱼渔业的可持续发展提供可供参考的科学依据。

通过气候变化表征因子相关性分析，得到了包含6个因子的相对独立的气候变化表征因子集{NAO, PDO、NPI、Niño1+2、dT、SOI}, 保证了神经网络模型输入变量的相对独立性。

通过因子互相关分析，确定了各气候变化表征因子与智利竹筴鱼渔获量相关性最强的滞后阶数，因子NAO、PDO、NPI、Niño1+2、dT以及SOI的最佳滞后阶数分别为4年、3年、3年、9年、12年与2年。

将6个气候变化表征因子的滞后数据以及智利竹筴鱼渔获量数据输入BP神经网络模型，以效率系数为评价规则得到了智利竹筴鱼预测最优模型。最优模型拟合效果显示，模型的渔获量拟合值与观测值有基本一致的变化趋势，线性相关系数为0.75, 拟合效果良好。

神经网络模型因子敏感性分析表明，在研究期间，影响东南太平洋智利竹筴鱼渔获量的气候变化表征因子主要为NAO、PDO和NPI。