秋刀鱼(Cololabis saira)为北太平洋中上层鱼类，是我国重要的大洋性经济种类之一，也是北太平洋渔业管理委员会优先管理鱼种之一^[1]。由于秋刀鱼的高度洄游性，在远洋渔业生产过程中，渔船每天需变更生产海域，并投入大量成本寻找渔场。但在目前国内外学者对秋刀鱼渔情预报的研究中，预报的时间分辨率多为月或年^[2,3,4,5], 针对周或天的渔场速报研究非常少。作为一种2龄的短生命周期鱼类，秋刀鱼的资源丰度和渔场分布易受到不同海洋环境的影响^[6,7], 因此，在预报模型建立中，各个环境变量的权重分析十分重要。

提升回归树方法可以有效计算出环境变量对模型建立的贡献大小^[8]。广义可加模型是多元线性回归的扩展，能够直接处理响应变量与多个解释变量的非线性关系^[9]。栖息地适宜性指数模型是评价环境对生物影响的重要方法，在渔业领域得到广泛的应用^[10], 该模型的兼容性强，能够融合BRT和GAM方法，既可以有效的选择环境变量，也能够很好的分析渔场与海洋环境间的复杂关系。本研究根据西北太平洋秋刀鱼生产数据和海洋环境数据，结合GAM和环境变量的权重分析，按月份构建秋刀鱼HSI模型，并分析模型在秋刀鱼渔场预报中的可行性，以期为秋刀鱼生产渔场速报提供参考。

渔业数据来自中国远洋渔业协会秋刀鱼技术组，渔船主要作业海域为38°~46°N, 150°~165°E, 数据时间为2013─2016年的7─11月，时间分辨率为天，数据项包括作业日期、经度、纬度、日产量(t)和作业船数等信息，数据总量为5620条。CPUE定义为每艘船每天的渔获量，单位为t/d, 作为渔业资源丰度的指标。采用渔场分布重心作为资源量的空间分布和变化指标^[11], 计算公式为：

式中，X、Y分别为渔场作业重心的经度和纬度，C_i为作业点i的CPUE, X_i和Y_i分别为作业点i的纬度和经度，n为作业总次数。

环境数据首先选用海表面温度(sea surface temperature, SST)及其变化海表温度梯度(sea surface temperature gradient, SSTG), 其通常作为判断渔场的关键指标^[6], 来源于美国国家海洋和大气管理局网站(http://www.noaa.gov), 其空间分辨率为0.01°×0.01°。SSTG利用梯度幅值公式求得^[12], 假设当前SST值为SST_i,j, 其相邻的4个网格点的SST分别为SST_i+1,j、SST_i-1,j、SST_i,j+1、SST_i,j-1,, 则

式中，i、j为正整数，分别表示网格数据行号和列号，$\Delta x$表示第j–1列与第j+1列之间经度方向的距离，$\Delta y$表示第i–1行与第i+1行之间纬度方向的距离，单位均为km。

其次是叶绿素a浓度(chlorophyll-a concentration, Chl-a), 其可以反映渔场的初级生产力条件，来自哥白尼海洋环境监测中心(http://marine. copernicus.eu), 空间分辨率为0.25°×0.25°。海面高度异常(sea level anomaly, SLA)和涡动能(eddy kinetic energy, EKE) 则可以作为动力环境和海流活跃程度的指示因子^[13], 来自卫星高度计资料，取自法国空间局AVISO网站(http://www.aviso. oceanobs.com), 其中包括海面高度数据和地转流数据，空间分辨率0.25°×0.25°, 涡动能的计算公式为^[13]:

式中，EKE为涡动能，单位为cm²/s²；${\rm{\;}}u' = u - \bar u$, $v' = v - \bar v$, $\bar u$和$\bar v$分别是纬向速度u和经向速度v的年平均值。最后，混合层深度(mixed layer depth, MLD)可以作为海水垂直分层的指示因子^[6], 来自全球海洋HYCOM模式数据(https://www.hycom. org), 空间分辨率1°/12×1°/12。本研究采用的环境数据时间分辨率均为天。

根据渔船作业的日期和经纬度信息，查找相应日期海洋环境数据网格上与该作业位置距离最近的经纬度所对应的数据值，匹配为该作业位置的环境数据经纬度和环境数据值。

BRT方法是在运算过程中随机抽取一定量的数据分析自变量对因变量的影响程度，而剩余数据用来对拟合结果进行检验，对生成的多重回归树取平均值输出，从而得出自变量对因变量的影响大小^[8]。本文BRT模型的数据计算采用了R软件中的梯度提升函数(gradient boosting machine, GBM)^[14], 其中抽样率(train. fraction)设置为0.8, 重复循环计算1000次，获取环境变量的权重分布，并取平均值计算标准差进行检验。

本研究利用GAM模型的样条平滑函数建立各个环境变量与SI的关系模型，表达形式如下：

式中，SI为变量的适宜性指数；$\alpha $为模型的截距；f为非参数平滑函数；${\rm{\;}}{x_i}$为第i\;个解释变量；$\;\varepsilon $为残差，$\varepsilon = {\sigma ^2}$且$E(\varepsilon ) = 0$。模型采用的是样条平滑法，误差分布估计为高斯分布^[15]。通常研究以CPUE建立SI, 假定最高CPUE为资源丰度最大的海域。但秋刀鱼的生产数据每天的差距量级可达100, 直接使用CPUE会降低SI指数精度。因此本文对每个环境变量进行数据分组，计算每个环境变量在各组距范围内的CPUE累积频率，建立每个组距下的SI, 计算公式为：

式中，SI_m为环境变量在每个组距范围内的适宜性指数，CPUE_m为每个组距范围下的CPUE累积频率，CPUE_max为每个组距范围下的最大CPUE累积频率。根据每个组距下的SI和环境数据进行样条平滑函数拟合，依据拟合优度${R^2}$, 确定各个环境变量的最适组距，并进行显著性F检验，当P<0.01时，则通过显著性检验，并依次分别求出各变量的SI指数。

结合权重分析和SI指数，本文利用赋予权重的算数平均算法建立HSI模型，算数平均算法具有不受SI极值影响的特点，是较为常见的有效方法，计算公式如下：

式中，i为第i个环境变量，n为环境变量总数，${W_i}$为第i个环境变量的权重，SI_mi为第i个环境变量的SI值。HSI取值范围一般为0~1, 0表示不适宜栖息地，1表示最适宜栖息地^[10]。本文利用2013—2015年7—11月的数据分月建立HSI模型，统计以0.1为间距的HSI值中产量占当月总产量的比重和平均CPUE, 对模型进行检验。然后将2016年7—11月的环境数据分别输入模型，得到HSI预报值，并绘制2016年生产作业点与HSI指数的空间分布图，以对渔场预报效果进行评价。

根据BRT方法分别计算各月环境变量SST、SSTG、Chl-a、EKE、SLA和MLD对CPUE的影响，结果显示，各个环境变量的权重标准差均小于0.05, 故可以使用平均值代表权重大小(表1), 不同环境变量对CPUE的拟合权重贡献率分布在10.9%~20.6%之间，7月最高和最低权重相差最大为9.4%。不同的月份权重分布各不同，8—10月排在前3位的变量分别为SSTG、SST、MLD, 7月为SSTG、MLD、EKE, 11月为SSTG、Chl-a、SST, 其中SSTG在秋季9—11月的权重值均为最高，而第4位之后各月权重排名均不相同。考虑到环境变量之间的相互作用会对模型产生影响，因此本文选择各月权重排在前4位的变量建立模型，以减少这种相互作用，并重新计算权重贡献率，结果见表2。

表1

表1 　筛选变量的权重 Tab. 1 Weight of selected parameter % 月份month W sst W sstg W chla W sla W eke W mld 7 14.25 19.20 16.35 13.45 16.46 20.30 8 16.99 19.76 10.88 15.95 16.13 20.30 9 17.88 19.30 15.87 14.12 13.62 19.21 10 18.36 20.56 14.16 15.73 14.99 16.19 11 17.18 17.87 17.61 15.57 14.67 17.10

表1 　筛选变量的权重 Tab. 1 Weight of selected parameter %

表2

表2 　各变量的SIm拟合结果 Tab. 2 Result of fitted SIm 月份month 环境变量parameter F R 2 组距步长step interval 权重/%weight 7 July SSTG 23.7 0.86 0.007 ℃/km 27.06 EKE 35.11 0.93 11 cm2/s2 22.91 Chla 10.21 0.69 0.08 mg/m3 22.71 MLD 22.22 0.68 2 m 27.32 8 August SST 117 0.98 0.6 ℃ 23.98 SSTG 42.22 0.81 0.002 ℃/km 26.16 EKE 104.5 0.97 11 cm2/s2 22.36 MLD 36.72 0.92 3.1 m 27.50 9 September SST 278.4 0.98 0.8 ℃ 24.88 SSTG 39.62 0.94 0.015 ℃/km 26.91 Chla 17.31 0.76 0.07 mg/m3 21.36 MLD 19.17 0.93 4 m 26.85 10 October SST 100.6 0.98 1 ℃ 25.21 SSTG 27.4 0.93 0.018 ℃/km 27.62 SLA 67.88 0.96 2.5 cm 24.41 MLD 52.86 0.97 5.9 m 22.76 11 November SST 56.83 0.98 1 ℃ 25.16 SSTG 181.9 0.98 0.022 ℃/km 25.07 Chla 30.65 0.77 0.15 mg/m3 25.06 MLD 17.89 0.90 5.8 m 24.72

表2 　各变量的SIm拟合结果 Tab. 2 Result of fitted SIm

结合权重分析结果，依据样条平滑函数拟合优度${R^2}$、自由度、F值以及显著性检验，循环计算得出7─11月各个环境变量的组距步长，结果见表2。结果显示，环境变量均在P<0.01水平下影响显著，通过显著性F检验。环境变量拟合方程的拟合度${R^2}$对于SST、SSTG、EKE、SLA均达到0.8以上，其中Chl-a在7月、9月和11月分别为0.69、0.76、0.77, 7月的MLD为0.68, 较高的${R^2}$说明GAM能够较好地拟合SI与各个环境变量的关系。

结合7─11月环境变量的SI拟合指数和各自对应的权重大小，运用算数平均值方法建立HSI模型，分别如下：

${\rm{HS}}{{\rm{I}}_7} = {\rm{S}}{{\rm{I}}_{{\rm{mld}}}} \cdot 27.32{\rm{\% }} + {\rm{S}}{{\rm{I}}_{{\rm{sstg}}}} \cdot 27.06{\rm{\% }} + $

${\rm{S}}{{\rm{I}}_{{\rm{eke}}}} \cdot 22.91{\rm{\% }} + {\rm{S}}{{\rm{I}}_{{\rm{chla}}}} \cdot 22.71{\rm{\% }}$

${\rm{HS}}{{\rm{I}}_8} = {\rm{S}}{{\rm{I}}_{{\rm{mld}}}} \cdot 27.5{\rm{\% }} + {\rm{S}}{{\rm{I}}_{{\rm{sstg}}}} \cdot 26.16{\rm{\% }} + $

${\rm{S}}{{\rm{I}}_{{\rm{sst}}}} \cdot 23.98{\rm{\% }} + {\rm{S}}{{\rm{I}}_{{\rm{eke}}}} \cdot 22.36{\rm{\% }}$

${\rm{HS}}{{\rm{I}}_9} = {\rm{S}}{{\rm{I}}_{{\rm{sstg}}}} \cdot 26.91{\rm{\% }} + {\rm{S}}{{\rm{I}}_{{\rm{mld}}}} \cdot 26.85{\rm{\% }} + $

${\rm{S}}{{\rm{I}}_{{\rm{sst}}}} \cdot 24.88{\rm{\% }} + {\rm{S}}{{\rm{I}}_{{\rm{chla}}}} \cdot 21.36{\rm{\% }}$

${\rm{HS}}{{\rm{I}}_{10}} = {\rm{S}}{{\rm{I}}_{{\rm{sstg}}}} \cdot 27.62{\rm{\% }} + {\rm{S}}{{\rm{I}}_{{\rm{sst}}}} \cdot 25.21{\rm{\% }} + $

${\rm{S}}{{\rm{I}}_{{\rm{sla}}}} \cdot 24.41{\rm{\% }} + {\rm{S}}{{\rm{I}}_{{\rm{mld}}}} \cdot 22.76{\rm{\% }}$

${\rm{HS}}{{\rm{I}}_{11}} = {\rm{S}}{{\rm{I}}_{{\rm{sst}}}} \cdot 25.16{\rm{\% }} + {\rm{S}}{{\rm{I}}_{{\rm{sstg}}}} \cdot 25.07{\rm{\% }} + $

$\;\;{\rm{S}}{{\rm{I}}_{{\rm{chla}}}} \cdot 25.06{\rm{\% }} + {\rm{S}}{{\rm{I}}_{{\rm{mld}}}} \cdot 24.72{\rm{\% }}$

通过对不同HSI组距下，秋刀鱼的产量比重和平均CPUE的计算，图1结果显示，2013─2015年7─11月的HSI值在0.5~0.9范围内的产量比重和分别为68.5%、75.6%、90.8%、86.6%、88.7%, 平均为82.0%, 其中秋季9─11月平均为87.7%；HSI值在0.5、0.6、0.7、0.7、0.7时产量比重最高，分别为25.5%、28.8%、27.7%、23.5%、25.5%, 各月最高CPUE对应的HSI值分别为0.8、0.6、0.6、0.7、0.7。由2016年7─11月的生产数据对模型进行评价，计算得到HSI值在0.5~0.9范围内的产量比重和，分别为99.4%、32.9%、67.3%、89.4%、76.9%, 平均为73.2%, 其中秋季9─11月平均为77.9%；HSI值在0.6、0.4、0.5、0.8、0.5时的最高，产量比重分别为27.5%、51.5%、35.0%、28.2%、34.1%, 各月最高CPUE对应的HSI值分别为0.6、0.4、0.5、0.8、0.4。总体来看，7─9月的预报结果与模型检验相差较大，8月的预报效果最差。秋季9─11月比夏季7─8月的HSI模型效果好，秋季为秋刀鱼的盛渔期，2013─2016年10月的平均CPUE达23.7 t/d, 而8月的平均CPUE为7.3 t/d。在秋季，以10月效果最佳，其预报HSI超过0.5以上时产量比重达89.4%, 且HSI值在0.2~0.8范围时，随着HSI指数的增加，CPUE呈逐渐上升趋势。本文以HSI值大于0.5为适宜秋刀鱼渔场形成海域，HSI值大于0.7为高值区，由以上结果可见，秋刀鱼高产渔区基本集中在HSI指数大于0.5的范围内，模型效果良好。

图1

图1 　不同HSI值下的产量比重和平均CPUE Fig. 1 Mean CPUE and percentages of yield under different HSI values

本研究利用各月的预报模型分别绘制2016年7─11月各月20─22日的HSI分布图，并叠加匹配各月每天的生产作业位置，结果如图2所示。7月HSI高值区主要分布在45°N以北和160°E以东海域，尽管HSI值平均在0.6左右，但实际生产点基本落在HSI高值范围，故7月预测的产量比重为最高月份。8月HSI高值区向西南方向移动，同时高值范围的面积逐渐减少，而渔船则向南移动越过了HSI高值区，作业位置完全分布于HSI小于0.5区域内，故8月预测的产量比重为最低月份。9月HSI高值区继续向南移动，高值面积集中且与渔船分布均呈条带状分布，作业位置集中落在高值区域，因此，尽管9月HSI值平均为0.5, 但高值范围内的产量比重高于8月。10月为秋刀鱼盛渔期，整个海域HSI高值海域面积范围广，作业位置集中在42°N和153°E附近的海域，且基本分布于HSI值大于0.8的范围内，因此预测结果最佳。11月HSI高值区主要分布于42°N和155°E以东的海域，而作业位置仍位于155°E以西的高值区域，故产量比重较10月低。总体上看，HSI高值区域和实际渔场的位置较为接近，且二者的空间分布移动方向基本一致。

图2

图2 　2016年7─11月各月20─22日秋刀鱼HSI和作业位置(红色圆圈)分布 Fig. 2 The distribution of HSI and CPUE in the fishing ground (red circle) of Pacific saury from July to November

作为短周期高度洄游鱼类，秋刀鱼对温度具有高度的敏感性，且其洄游路线很长。一般在夏末秋初，秋刀鱼从亲潮水域南下向黑潮水域移动，幼鱼逐渐成熟，秋季移动至黑潮附近，黑潮暖水与亲潮冷水交汇，海洋动力过程活跃，往往形成温度锋面现象，引起海水垂直混合加强，营养盐富集并伴有浮游植物大量繁殖，形成高生产力区，从而形成秋刀鱼盛渔期渔场^[16]。在整个洄游过程中，秋刀鱼对SST的响应更多呈现出季节性变化的特点。而在盛渔期，鱼群集聚，捕捞作业范围缩小，以往研究也发现，秋刀鱼作业海域的温度锋面越活跃，作业位置距离锋面越近，单位渔获量越高^[17,18], 因此，在秋季SSTG对渔场分布的影响更加突出。本研究以时间分辨率为月份的统计分析中，SSTG和SST的权重值排序与秋刀鱼的洄游集群特点基本一致。

其次，本研究发现，MLD的权重贡献率基本在前3位，这与秋刀鱼对SSTG的响应特征情况相似。由于秋刀鱼的垂直洄游特性，即白天分布水深为30~40 m水层，夜间为0~5 m的近表层^[19], 而在温度锋面处由于海水的垂直混合从而MLD加大，这一变化对应着秋刀鱼可适应的上层水体向下延伸，从而扩展秋刀鱼的上下活动范围。因此，MLD的变化影响鱼群的栖息水层分布，进而影响了鱼群集聚和渔场分布。第三，对于Chl-a、SLA和EKE 3个环境变量，主要体现在叶绿素、海流、涡旋等对鱼群的影响，在不同的月份由于鱼群成熟期不同，因而影响各不相同；在短期尤其盛渔期期间，可能与这3个变量在渔场区域的变化幅度相对于SST、SSTG、MLD较小有关，因而对渔场分布的影响权重贡献率略低。

本研究利用GAM方法的样条平滑函数结合环境变量的权重分析结果，分月建立了HSI模型，并分别对模型进行了检验和评价，总体精度分别为82.0%和73.2%, 其中秋季可达87.7%和77.9%。但7月和8月模型检验和评价结果相差较大，导致这个现象的原因可能是由于，在夏季，秋刀鱼处于成长期，在洄游过程不断改变方向寻找适宜温度等适宜环境，中心渔场不确定，这导致渔船需要花费大量时间寻找渔场，当遇到捕捞产量减少时，出于经济效益考虑，渔船往往随之迅速转移渔场，使得实际作业时间较短，渔场位置也在不断变化中。由于本文采用生产数据，作业位置一般会聚集在较小的空间区域，可能导致局部环境变量的权重被放大，从而降低了非传统作业区域的预报精度。此外，不同年份下的大尺度气候模态不同，气候模态通常直接影响温度的位置变化，会导致渔场在这一年整体偏南或者偏北^[5], 因此，本研究的模型在夏季预报中不稳定，对于7月和8月的预测并不适用。但模型在秋季的预测效果良好，尤其是10月，预测精度达89.4%, 这是因为在盛渔期，秋刀鱼在温度季节性变化以及其洄游特性的影响下，向南迁移集中至黑潮区域，这一区域环境变量(SSTG、MLD等)的改变有利于秋刀鱼的生存，本研究中，权重分析的结果亦是反映出了这一变化结果，因此，该模型在短期速报中具有优势，适用于秋季集中作业期间，能更快地确定渔场方向和范围。

相较于以往神经网络、最大熵以及基于灰色系统等西北太平洋秋刀鱼预测模型来说^[6,20,21], 本研究通过HSI模型的算法很好的融合了BRT和GAM方法，体现出了不同环境变量的重要性以及减少了各个环境变量之间的相互影响，且该方法灵活易操作，具有较强的实用性。但本研究的栖息地预报准确度总体而言仍不够精确，模型在盛渔期每天的预测中具有优势，在其他时期有一定的局限性，需要进一步将大尺度环境变量如厄尔尼诺、南方涛动等加入分析，以及根据秋刀鱼的洄游和集群规律，进一步开展不同时间的局部海洋动力环境特征研究，以提高西北太平洋秋刀鱼渔场速报精度，提高远洋渔业生产效率和经济效益。