山东南部近海秋季长蛇鲻相对资源量的分布及其影响因素❋

刘逸文， 张崇良❋❋， 昝肖肖， 孙 铭， 徐宾铎， 薛 莹， 任一平,3

(1. 中国海洋大学水产学院, 山东 青岛 266003； 2. 山东大学(威海)海洋学院， 山东 威海 264209；3.青岛海洋科学与技术试点国家实验室 海洋渔业科学与食物产出过程功能实验室, 山东 青岛 266237)

长蛇鲻(Sauridaelongate)属暖温性近海底层鱼类，在中国渤海至南海均有分布，是中国重要的经济鱼种之一。长蛇鲻的相关研究主要集中于年龄、生长、食性、遗传学等方面[1-3]，有关其资源分布的研究较少，影响其分布的关键因素尚不明确。许多研究表明，在个体发育的不同阶段，鱼类的分布模式会发生变化[1, 4]，成体和幼体的分布与影响因子的关系可能存在差异。因此探究成体与幼体的环境适应性差异，有助于解析鱼类的分布模式及变化，识别和保护重要的鱼类栖息地[5]，对于渔业资源的合理利用与保护有着重要意义。

当前物种分布的研究越来越多地依赖于统计学方法，利用统计模型评估物种分布与影响因子的关系，其中，广义可加模型(Generalized additive model，GAM)是研究物种分布常用的一种统计学模型。广义可加模型是广义线性模型的非参数化扩展，其优点是能够处理高维数据中响应变量与解释变量间的非线性、非单调关系[6]。在渔业资源研究应用中，水温、盐度、水深等自变量对于资源密度等因变量的影响往往不是线性的，因此，GAM模型的应用有十分重要的意义[7]，已被越来越广泛地应用于渔业资源研究[5, 7-10]。

本文根据2016—2017年在山东南部海域进行的渔业资源底拖网调查数据，利用GAM方法重点分析了秋季长蛇鲻的资源分布及其影响因素，并探讨了影响长蛇鲻幼体与成体分布差异的原因，旨在更深入的了解长蛇鲻幼体和成体的生态习性和分布规律，揭示其生活史差异导致的环境适应性差异，为长蛇鲻资源的合理利用和养护提供科学依据。

1 材料与方法

1.1 数据来源

研究数据来源于2016年10月、2017年1、5、8月在山东南部进行的4个航次(见表1)的渔业资源调查，调查范围为119°30′E～124°E、35°0′N～36°45′N，共设置63个站位(见图1)。用于样品采集的调查船为单拖底拖网渔船，主机功率为220 kW，平均拖速为3 kn，拖网时间1 h，拖网时网口高度约7.53 m，宽约15 m。将渔获物冷冻保存，带回实验室内进行样品鉴定和生物学测量。调查及分析均按《海洋生物调查规范》(GB 12763.6—2007)[11]要求进行。

图1 山东南部近海渔业资源底拖网调查站位

调查发现，本研究海域长蛇鲻由于繁殖和洄游[1]等的影响，分布具有明显的季节性变化，在春、夏和冬季航次捕获率很低(见表1)。长蛇鲻的繁殖期为5～7月[1]，春、夏航次均无长蛇鲻幼体的分布，而夏季长蛇鲻虽资源量较高但呈现集群分布，仅在小部分站位出现。冬季长蛇鲻群体进行越冬洄游，于调查海域分布较少，相对资源量和出现频率均为最低。秋季长蛇鲻的相对资源量和出现频率均为年内最高，满足模型的数据量要求，因此本文仅选用2016年10月(秋)数据进行分析。

本研究中涉及的环境因子包括水温、盐度、水深、叶绿素a含量、离岸距离和底质类型。其中，拖网时同步使用CTD仪(CTD75M/1167)测量各个站位的水温、盐度、水深、叶绿素a含量等环境因子。离岸距离根据站位经纬度计算得来，底质类型包括砂、砂-粉砂-黏土、淤泥质粉砂、粉砂质黏土四类。

本研究同时考虑生物因子(饵料生物丰度)[12]对长蛇鲻分布的影响。生物因子的选择基于已发表的长蛇鲻摄食生态研究[2]，包括戴氏赤虾(Metapenaeopsisdalei)、短鳍鱼衔(Callionymussagitta)、六丝钝尾虾虎鱼(Amblychaeturichthyshexanema)、枪乌贼(Loligospp.)、鳀(Engraulisjaponicus)和细条天竺鲷(Apogonlineatus)。

表1 山东南部近海长蛇鲻的相对资源量

注：平均值±标准误差。Note: Mean ± standard error.

1.2 数据处理

将各站位长蛇鲻及饵料生物的实际渔获量标准化为拖速3.0 kn与拖网时间1.0 h的相对资源量(g/h)，作为资源量指标。根据本次调查的生物学测量数据，以全长197.5 mm作为50%性成熟体长，并根据50%性成熟体长将长蛇鲻划分为幼体和成体[13]。

在构建模型前对环境因子进行Pearson相关性分析，以剔除相关性较高的因子对模型拟合的影响。结果显示，表温和底温，表盐和底盐，水深、底温和离岸距离间存在显著的相关性，可能对模型拟合造成影响。根据长蛇鲻底栖的生态习性，初步筛选出5个因子：水深、底层盐度、底层水温、表层叶绿素A含量、底质类型。

根据长蛇鲻与饵料生物的空间重叠指数(Schoener indices, SIs)[14]筛选生物因子，公式如下：

(1)

式中：Zi(x)和Zi(y)分别表示x、y物种在i站位的相对资源量；n为站位个数。结果如表2所示，仅有枪乌贼的重叠指数较高，为0.41，其它饵料生物的重叠指数范围为0.01～0.18。本文剔除重叠指数小于0.1的物种，筛选出枪乌贼、细条天竺鲷、短鳍和六丝钝尾虾虎鱼作为生物因子，将其相对资源量进行对数化处理作为丰度指标加入到模型中。

表2 饵料生物与长蛇鲻的空间重叠指数

1.3 GAM模型

以长蛇鲻相对资源量作为响应变量，利用GAM模型分析响应变量与解释变量(环境因子和生物因子)之间的关系。GAM模型的表达式如下[6]：

(2)

式中：g(μ)为连接因变量的连接函数；xi为解释变量；α表示模型中的截距；ε为随机误差项，fi(xi)为解释变量的非参数平滑函数，通过样条平滑得到。

采用向前逐步筛选的方法构建模型，首先将初步筛选的9个因子逐个分别带入模型，根据AIC(Akaike information criterion)准则[15]，选择AIC最小的单因子模型，逐渐加入其他因子，筛选出AIC值最小的双因子模型，循序渐进直至得到一个最小AIC的模型，即为拟合效果最好的模型[16]。

AIC的计算方法如下[15]：

(3)

模型筛选过程中，考虑到水深与底温之间显著的相关性，对底温进行残差处理[17]，以排除解释变量间的共线性影响：

res=residual[lm(sbt～depth)] 。

(4)

式中：sbt代表底温；depth代表水深；lm表示线性回归。在底温和水深同时存在的模型中，以res底温残差值代替底温进行模型的拟合。

本研究对成体、幼体的分布分别进行了建模和检验，其方法是首先拟合了分类模型M1，分别模拟成体和幼体相对资源量与影响因子的关系。其中平滑函数拟合了对成体、幼体不同的回归曲线，即与影响因子的不同关系。此外以相同数据集拟合集合模型M2，在同一站位中成体、幼体相对资源量相加进行模型拟合，即反映了总资源量与影响因子的共同关系。

利用ANOVA方差分析，对M1、M2两个嵌套模型的差异进行检验，根据卡方检验结果验证两者间是否存在显著差异，以检验成体和幼体与影响因子关系的差异显著性。

1.4 模型检验

利用交叉验证(Cross validation)检验所得模型的预测效果，选取80%作为训练数据构建模型，剩余的20%作为验证数据，与模型的预测结果相对照。交叉验证过程运算1 000次，通过计算预测值与实际值间的均方根误差(Root Mean Squared Error，RMSE)、线性关系截距斜率和决定系数R2等4个参数判断模型的预测效果。RMSE用来衡量模拟值和真值之间的偏差，决定系数R2表示可由预测变量变异解释的预测值线性变异的比例，截距描述模型的偏差，表示观测与模拟之间系统恒定的差异，斜率描述模型的一致性，表示观测值和实际值之间的差异，解释为低估或高估[18-19]。

1.5 分布预测

本研究使用两阶段的GAM模型[12]进行长蛇鲻资源分布的预测，一阶段仅根据环境因子构建饵料生物(Xj)的分布模型：

(5)

将一阶段的结果，即饵料生物丰度作为变量加入二阶段GAM模型，预测长蛇鲻的相对资源量(Y)分布：

(6)

以上两式中：xj为环境因子；Xj为饵料生物j的丰度；Y为长蛇鲻的相对资源量。

用于预测的环境数据提取自FVCOM(Finite-Volume Community Ocean Model)模型，FVCOM是一种无结构网格、有限体积、自由表面三维原始方程的海洋环流模型[20]。本文根据山东南部海域的FVCOM模型，提取33 825个数据点的环境数据(水温、盐度、水深)预测长蛇鲻的资源分布，并绘制相对资源量分布图。

利用Surfer13软件绘制长蛇鲻相对资源量分布图，利用R(version 3.4.4)中的mgcv包[21]完成模型的拟合和检验过程。

2 结果与分析

2.1 长蛇鲻资源分布特征

长蛇鲻相对资源量的分布在成体、幼体间存在明显差异。秋季幼体的分布范围为119°30′E～122°15′E、35°N～36°30′N，主要分布在30 m等深线附近及以浅水域；成体分布范围为119°30′E～122°30′E、35°N～36°45′N，较幼体分散，分布水深范围更广，大致存在海州湾水域及调查水域中部两个分布中心(见图2)。

(左：2016年10月成体；右：2016年10月幼体。 Left: Adult in October 2016； Right: Juvenile in October 2016.)

2.2 GAM模型

根据AIC准则筛选出的最优模型为：

log(Y+1)=α+s(log(lsp+1))+s(res)+s(depth)+s(sbs)+ε。

(7)

式中：Y为长蛇鲻的相对资源量，对Y进行对数化处理得log(lsp+1)作为响应变量；lsp为枪乌贼相对资源量；res为残差处理后的底层水温；depth为水深；sbs为底层盐度。区分和不区分成体与幼体的两模型M1和M2对应的AIC值分别为523.17、545.07，累计偏差解释率分别为78.3%、69.3%(见表3)。方差分析结果表明，M1、M2两模型差异极显著(P<0.01)，即成体、幼体与影响因子的关系有极显著差异。

方差分析表明，对长蛇鲻成体相对资源量有显著影响的因子包括枪乌贼丰度、底层水温和水深，对长蛇鲻幼体相对资源量有显著影响的因子包括乌贼丰度、底层盐度和水深，对长蛇鲻整体资源量影响显著的因子为底层水温和水深和底层盐度(见表3)。

表3 GAM拟合结果的方差分析

Note: lsp = relative abundance ofLoligospp., sbt = sea bottom temperature, sbs = sea bottom salinity.

GAM模型表明(见图3)，长蛇鲻成体的相对资源量随枪乌贼丰度和底层水温的增加均呈现先上升后下降的趋势，幼体相对资源量随枪乌贼丰度和底层水温的增加均呈现上升趋势；成体、幼体的相对资源量随水深增加的变化有所差异，成体在25～35 m的区间内相对资源量急剧下降，随后缓慢上升进而下降，在45～50 m水域附近出现小的峰值，幼体相对资源量随水深的增加始终保持稳定下降的趋势；成体、幼体的相对资源量随盐度的变化趋势明显不同，成体的变化趋势不明显，幼体的相对资源量随盐度的增加呈现稳定上升的趋势。长蛇鲻整体的相对资源量随枪乌贼丰度和底层水温的增加均呈现先上升后下降的趋势；随着水深的增加，相对资源量呈现总体下降的趋势，其中，在35～50 m出现小的峰值；相对资源量随盐度的增加呈现上升趋势。

图3 M1和M2中解释变量对长蛇鲻相对资源量的影响

2.3 长蛇鲻相对资源量的预测

交叉验证的结果表明(见表4)，M2模型有较低的均方根误差和较高的决定系数，两模型线性回归平均的截距和斜率差异不大，M2更近于无偏估计，且预测值的波动范围更小(见图4)。因此，不区分成体幼体的M2模型具有更好的预测效果。

表4 模型交叉验证结果

(虚线表示1 000次线性回归的平均值，实线代表模拟值与真值的无偏估计。The dotted lines represent the mean of 1000 times linear regressions. The solid lines represent the unbiased prediction between simulated values and true values.)

图4 模型1 000次交叉验证的散点图和拟合曲线

Fig.4 Scatter plot and fitting curves for results of 1 000 times cross validation

基于预测的枪乌贼分布和提取自FVCOM的环境因子数据，本文根据已构建的M1模型，绘制了山东南部近海秋季长蛇鲻成体和幼体的资源分布图(见图5)。长蛇鲻的分布模式与实际观测接近，成体存在海州湾和调查海域中部两个分布中心，在深水区的分布较少；幼体的分布呈现明显的随水深分层的现象，主要分布于30 m等深线以浅的区域，近岸相对资源量高，在深水区无分布。

图5 山东南部近海秋季长蛇鲻成体和幼体的相对资源量分布图

3 讨论

许多鱼类的栖息分布随着生活史不同阶段发生变化[22-23]，而本文作为一个研究案例，揭示了长蛇鲻成体和幼体不同的分布模式。相对资源量分布图中长蛇鲻成体主要分布于60 m以浅海域，在123°E以东的海域几乎无分布(见图5)。可能由于运动能力发达，长蛇鲻成体有着广泛的分布，在海州湾及整个沿岸带形成了第一个高值区，在调查海域中部形成了第二个高值区。幼体分布范围相对狭小，呈现近岸高远岸低的主要趋势(见图5)，这也表明近岸海域是长蛇鲻的重要的产卵场和育幼场，这与以往的研究结果相一致[24]。

本文分析了生活史差异对模型拟合的影响，两模型存在极显著差异(P<0.01)，表明成体与幼体的分布随影响因子的变化趋势存在极显著差异。不区分成体、幼体的模型M2会忽略成体、幼体间存在的部分差异(见图3)。因此，在模型拟合的过程中，需要考虑个体不同生长阶段的差异，以增加模型的准确性。本研究也发现，相比于M2模型来说，考虑到生活史差异的M1模型有较大的预测误差，这也说明随着模型复杂度的增加，模型的预测效果可能下降，预测误差可能增加[25]。因此，模型构建应考虑到拟合效果和预测效果间的权衡问题，避免过度拟合的出现。

长蛇鲻的相对资源量与生物因子和环境因子(水温、水深、盐度)显著相关。其中饵料生物通常表现为正效应，本研究中随枪乌贼丰度的增加，长蛇鲻成体和幼体的相对资源量均呈现上升的趋势，但在枪乌贼丰度达到e8(g/h)后，成体的相对资源量出现下降的趋势，这可能说明在较高饵料密度条件下其它因素限制了长蛇鲻成体的资源量。水温是影响鱼类生长、发育和分布的主要环境因子之一，有研究表明长蛇鲻的适温区间为20.9～28.4 ℃[26]。本研究中，海域底温范围为7.1～22.3 ℃，成体相对资源量随水温增加呈先上升后下降趋势，推测成体的最适水温低于22.3 ℃，幼体始终为上升趋势，推测幼体的最适水温高于22.3 ℃，即成体幼体的最适水温存在差异。Xue等的研究表明，底温在>21.8 ℃的区间内开始引起长蛇鲻相对资源量下降[12]，本研究结果与其结论相近。水深对长蛇鲻相对资源量表现为负效应，随水深的增加，成体与幼体相对资源量均表现为下降趋势。研究发现，在55 m以深水域无长蛇鲻的分布，这可能是由于调查季节受到南黄海东侧冷水团中心的影响，底层水温低于10 ℃[27]，不适合长蛇鲻栖息分布。由于本研究中水深与底层水温有较强的负相关性，这种分布模式可能是水温与水深共同作用的结果，对成体和幼体的分布均表现出限制作用。另外长蛇鲻相对资源量也受到盐度的影响，盐度的变化对成体相对资源量无显著影响，但对幼体影响极显著，随底层盐度的上升，幼体相对资源量有明显增加，这与Xue等的研究结果相似[12]。综上所述，本研究认为成体和幼体间的主要差异表现在对水温和盐度的适应性，幼体可能对盐度的变化更敏感，而成体受水温的影响可能较大。

与传统回归模型相比，GAM模型可以综合多个变量，定量描述资源与影响因素间的关系，评估各变量的重要程度，并可对资源密度等进行预测[28-29]，本研究中建模过程可以为相关研究提供有益经验。例如模型的回归方程存在多个解释变量时，变量间可能具有不同程度的共线性，强共线性会使模型的拟合中出现方差膨胀而增大误差，因此在模型构建前应加以排除[30-31]。本研究数据中，底层水温和水深间具有较强的相关性，考虑到水深是一个综合因子，为保留更多原始数据信息，本文将水深也加入分析，在模型拟合过程中，对底层水温进行残差处理，以res值代替水温进行拟合，去除相关性的影响。这种处理得到的结果可以理解为水深保持一定时，相对资源量随底层水温的相对变化趋势，但res不代表实际的数值，分析长蛇鲻的最适水温区间等信息存在一定困难。

GAM模型研究要求高准确度、可信度以及全面的调查数据，以尽可能包含影响物种分布的所有因素。本文考虑的环境因子包括水温、盐度、水深、叶绿素a含量、底质类型，但未能包含所有重要环境因子，如溶解氧等。此外，长蛇鲻作为一种底层鱼类，底质类型可能是影响其分布的原因之一，但在模型的拟合过程中，底质类型并没有表现出对相对资源量显著地影响，在最优模型中予以排除。这可能是由于调查海域的底质类型均一性较强[32]，对长蛇鲻的栖息分布影响不大。此外本文选定的底质类型分类方式较为粗糙，只根据物理特性将底质分为砂质、黏土类和砾石类，而忽略了其生物因素(如植被覆盖率等)。另一方面，调查海域是长蛇鲻秋季的重要栖息地，在特定时期长蛇鲻相对资源量大、分布广泛，是鱼类群落中重要的优势种[12]，确有必要进行物种分布研究。因此在今后的研究中，需进一步优化研究的时间、空间尺度，覆盖长蛇鲻的分布区域，在更大的空间尺度上研究长蛇鲻不同季节的分布模式，并获取更全面的环境数据，增加模型的准确性和预测效果。