基于GAM 的黄海南部越冬小黄鱼资源丰度与环境因子关系

李国东，李冬佳，熊 瑛，仲霞铭，宋大德，康中杰，汤建华，施金金，吴 磊

（1.江苏省海洋水产研究所，江苏南通 226007；2.上海海洋大学海洋科学学院，上海 201306）

小黄鱼（Larimichthyspolyactis）隶属鲈形目，石首鱼科，黄鱼属，系暖温性底层鱼类［1］，广泛分布于中国渤海、黄海、东海及朝鲜半岛西岸海域［2］。小黄鱼在我国渔业中占有重要地位，其资源曾经历兴盛期、衰退期和恢复期几个阶段［3-5］，海洋捕捞产量从20世纪80年代末期的2万t已上升到2020年的近30万t［6-7］。但与之相伴随的是小黄鱼低龄化、小型化和性成熟提前等特征日益凸显［3，8］。小黄鱼是一种洄游性鱼类，一年四季在产卵场、索饵场和越冬场之间有规律地做季节性洄游［9］。刘效舜等［10］研究表明，黄海存在3个主要小黄鱼种群：黄海北部-渤海种群、黄海中部种群以及黄海南部种群。春季黄海南部小黄鱼种群从沙外渔场迁徙到大沙渔场中西部，秋季则开始向东迁移到外海越冬场，在冬季抵达沙外渔场西部水域进行越冬［11］。鱼类越冬洄游是鱼类种群动态的重要组成部分，与海洋环境因子间有着密切联系。近年来，受气候变化和人类捕捞活动的影响，小黄鱼生长洄游习性发生了变化，越冬场位置有一定程度偏移［12-13］。因此，分析越冬小黄鱼资源丰度对环境变化的响应，对于研究小黄鱼种群补充和渔业资源管理具有重要意义。

核密度估计（kernel density estimation，KDE）是一种用于计算要素在其周围邻域中密度的非参数方法［14］。相较于其他空间分析方法，核密度估计法需要参数较少，受人为主观因素影响较小，因此是探究渔业资源空间分布模式的有效工具［15］。广义可加模型（generalized additivemodel，GAM）是一个能很好拟合高维数据中因变量与自变量之间非线性关系的生态位模型［16］，已被广泛应用于渔业资源分布与环境关系相关研究［17-18］。然而，关于越冬期间黄海南部小黄鱼种群资源分布与环境的关系尚少见研究报道。另外，目前近海渔业资源分布与环境因子关系的研究主要是利用调查站点的捕捞渔获数据及在该站位使用CTD等仪器测定的环境数据，如海水温度、海水盐度、叶绿素a浓度等［19-21］。随着遥感卫星技术的发展，海洋遥感环境数据精度不断提高。与使用CTD等仪器测量站点数据相比，遥感环境数据具有覆盖面广、人工成本低和自动化获取等诸多优势，促进了渔场分析预报等应用研究的发展［22］。

本文利用2010—2019年冬季黄海南部小黄鱼渔获数据，基于核密度估计法分析越冬期间小黄鱼资源分布情况，并结合相应时空的海洋卫星遥感监测数据，运用GAM探究黄海南部越冬小黄鱼资源丰度与环境因子的关系，以期为预测小黄鱼越冬渔场、制定小黄鱼资源合理捕捞和科学养护政策提供科学依据。

1 材料与方法

1.1 数据来源

本文小黄鱼捕捞数据来源于江苏省海洋水产研究所2010—2019年冬季（12月—2月）在黄海南部进行的帆张网商业渔船捕捞日志，主要包含捕捞小黄鱼所在渔区、小黄鱼产量以及所布网次数据。依据小黄鱼捕捞样本数据所在渔区（定义为10′×10′的空间分辨率）进行了地图可视化处理（图1）。

图1 2010—2019年捕捞越冬小黄鱼分布位置Fig.1 Distribution locations of overw intering Larim ichthys polyactis fishing data during 2010—2019

依据小黄鱼资源分布与环境因子关系相关研究［20-21］，初步选择水深（water depth，Depth）、底层温度（sea bottom temperature，SBT）、底层盐度（sea bottom salinity，SBS）和叶绿素a浓度（chlorophyll-a concentration，Chla）作为环境因子。其中，水深数据来源于英国海洋学数据中心（https：／／www.gebco.net／data＿and＿products／gridded＿bathymetry＿data），空间分辨率为36″，数据更新时间为2017年12月；叶绿素a浓度数据来源于NOAAOceanWatch（https：／／oceanwatch.pifsc.noaa.gov／），空间分辨率为3′，时间分辨率为月；底层温度和底层盐度数据来源于中科院大气物理所网站（http：／／www.ocean.iap.ac.cn），空间分辨率为30′，海水分层测定的垂直分辨率约为5～10 m，时间分辨率为月。

1.2 数据处理与分析

1.2.1 标准化处理

由于2010—2019年商业捕捞渔船数据及其携带的网具规格不一，需要对所有渔获量数据进行标准化处理，剔除网具规格差异导致的误差，因此采用单位捕捞努力量渔获量（catch per unit effort，CPUE）表示小黄鱼越冬资源丰度，计算方法如下：

式（1）中，C为小黄鱼渔获量（kg），E为该渔获量所投入的捕捞努力量，S为所使用的帆张网网具迎流面积（m2），N为投网次数。

由于捕捞数据和环境因子数据的空间分辨率不完全一致，且它们之间呈独立状态，需利用ArcGIS10.2软件进行处理。处理过程如下：通过捕捞渔获位置生成相应捕捞点集，基于此点集对各环境栅格数据采用最邻近分配法进行采样，并通过空间位置信息连接捕捞点集的资源丰度和各环境要素数据，最后进行点集合并及统计分析。

1.2.2 核密度估计法分析

本文基于计算得到的小黄鱼产量数据，运用核函数［14］进行分析，能直观反映出黄海南部越冬小黄鱼相对资源丰度的时空分布情况。

式（2）中，K［］是核函数；h为带宽；n是在带宽范围内的已知点数目；d是数据的维度。

1.2.3 影响因子筛选

本文选取小黄鱼越冬资源丰度作为因变量。除了初选的环境因子外，亦加入时间（年份）作为初始解释变量，探究黄海南部越冬小黄鱼资源丰度的年际变动情况。本研究采用条件数k［23］和方差膨胀因子（variance inflation factor，VIF）［24］对初始因子进行多重共线性检验，从而剔除存在共线性的因子。

设x1，x2，x3，…，xp是选取的初始因子，X1，X2，X3，…，Xp是经过中心化和标准化的向量，记X=［x1，x2，x3，…，xp］，设λ为XTX的一个特征值，φ为对应的特征向量，其长度为1，即φTφ=1。度量多重共线性严重程度的一个重要指标是矩阵XTX的条件数k：

式中，λmax（XTX）、λmin（XTX）分别表示矩阵的最大、最小特征值。条件数表征了XTX特征值差异的大小。若条件数k＜100，则表示初始因子之间存在较小的多重共线性程度；若100≤k≤1 000，则表示多重共线性为中等或略强程度；若k＞1 000，则表示多重共线性十分严重［23］。条件数k的计算可通过R语言kappa函数实现。

此外，VIF的平方根表示变量回归参数的置信区间能膨胀与模型无关的预测变量的程度，可通过使用R语言中car包中的VIF函数确定VIF值。一般认为VIF＞4，即存在多重共线问题［24］。

1.2.4 GAM模型的构建

本研究以小黄鱼资源丰度为响应变量，以初步筛选的相关因子作为解释变量，利用GAM模型对小黄鱼越冬资源丰度与相关环境因子的关系进行拟合。GAM模型的一般形式如下［25］：

式（4）中，Y为小黄鱼资源丰度，即每1 m2网具迎流面积的渔获量，单位为g·m-2；α是适合函数的截距；ε是服从正态分布的随机误差项；xj表示解释变量，即各站位的影响因子；fi（xj）是各环境变量xj的任意单变量函数，通过样条平滑函数来估计。模型的误差分布估计为Gaussian分布，连接函数为自然对数。

依照赤池信息准则（Akaike information criterion，AIC），在AIC值最小的因子预测函数上依次加入其他因子，进而得到AIC值最小的多因子预测模型［25］。其中，AIC值越低，模型的拟合效果越好；方差解释率越高，模型拟合效果越好。AIC的计算过程如下：

式中，k是参数个数，L是似然函数。

2 结果与分析

2.1 越冬期间小黄鱼资源分布

黄海南部越冬期间小黄鱼资源分布呈现一定的聚集性，主要分布在32°24′～34°00′N、123°00′～125°00′E区域。而各月份小黄鱼资源分布情况有所不同（图2～图4），其中12月小黄鱼资源核密度高值区域主要集中在32°48′～33°48′N、123°00′～124°10′E，即黄海南部水深45～60 m海域。1月小黄鱼资源逐渐迁移至水深55～75 m海域，集中在32°24′～34°00′N、123°24′～125°00′E区域，呈现东南-西北分布。2月小黄鱼资源分布区域与1月份基本相同，但资源核密度高值区域主要集中在大沙渔场中北部，且该区域的越冬小黄鱼资源核密度与1月相比明显增加。

图2 12月黄海南部越冬小黄鱼资源核密度分布Fig.2 Kernel density distribution of Larim ichthys polyactis resources in December

图3 1月黄海南部越冬小黄鱼资源核密度分布Fig.3 Kernel density distribution of Larim ichthys polyactis resources in January

图4 2月黄海南部越冬小黄鱼资源核密度分布Fig.4 Kernel density distribution of Larim ichthys polyactis resources in February

2.2 小黄鱼越冬洄游影响因子

通过条件数k和VIF对初始因子进行多重共线性检验（表1），结果显示初始因子的条件数小于100，VIF值均小于4，即初始因子的多重共线性程度均较小，因此构建GAM时将年份、底层温度、底层盐度、叶绿素a浓度以及水深作为本研究的解释变量。

表1 初始因子的多重共线性检验Tab.1 Multicollinearity test of influencing factors

通过GAM模型拟合小黄鱼越冬资源丰度与各影响因子间的关系，利用AIC准则进行因子筛选（表2）。根据AIC值最小原则得到最优GAM模型，模型的最终表达式为：

表2 最优GAM 模型中各影响因子的相关参数Tab.2 Relevant parameters of each factor in the optimal GAM model

式（6）中，Y为小黄鱼越冬资源丰度，α为截距，s（Year）为年份因子效应，s（SBT）为底层温度效应，s（SBS）为底层盐度效应，s（Chla）为叶绿素a浓度效应，s（Depth）为水深效应。

小黄鱼最优GAM模型及各影响因子相关参数表明，影响黄海南部小黄鱼资源越冬洄游的变量为水深、底层温度、底层盐度、叶绿素a浓度和年份。模型累积偏差解释率为54.4%。其中在各海洋环境因子中，对小黄鱼越冬洄游影响最大的因子为水深（34.3%），最小的为叶绿素a浓度（4.6%）。

2.3 小黄鱼越冬对环境因子的响应

本文建立GAM模型时加入年份因子作为初始解释变量，结果显示2010—2017年黄海南部越冬小黄鱼资源丰度的年际变动较小，而2017年后小黄鱼资源丰度有所上升（图5）。在选取的影响因子中，水深对小黄鱼的资源丰度影响明显最大（表2）。随着水深的增加，小黄鱼资源丰度呈现先减少后增加的趋势。其中，水深在30～40 m范围内，小黄鱼资源丰度随水深增加而减少；在40～70 m水深范围内，小黄鱼相对资源丰度总体呈现随着水深增加而明显增加趋势。底层温度对小黄鱼越冬资源丰度影响波动较大，但整体上资源丰度随着温度升高而振荡上升。底层盐度范围在32.8～33.5时，小黄鱼资源丰度随着盐度的增加而明显增加；当盐度大于33.5时，小黄鱼资源丰度一直较为稳定。叶绿素a浓度范围在1.0～2.3 mg·m-3时，小黄鱼资源丰度随着叶绿素a浓度增加而明显增加；当叶绿素a浓度大于2.3 mg·m-3时，小黄鱼资源丰度随着叶绿素a浓度增加而缓慢减少。

图5 基于GAM 模型的小黄鱼资源丰度对环境因子响应情况Fig.5 Response of Larim ichthys polyactis abundance to environmental factors based on GAM model

3 讨论

小黄鱼越冬洄游期间资源丰度高值区域一定程度上可看作为小黄鱼的越冬场。越冬场对于补充小黄鱼种群十分关键，因为越冬场汇集了一年中的大部分小黄鱼成鱼［4］。本研究表明，2010—2019年黄海南部越冬小黄鱼群体资源丰度高值区域在空间分布上均呈现一定的聚集性（图2～图4），主要分布在32°20′～34°00′N、123°00′～125°00′E区域，与之前研究结果基本相同［10-13］。此外，越冬期间各月份小黄鱼资源分布情况有所不同，其中12月小黄鱼资源核密度高值主要集中在32°48′～33°48′N、123°00′～124°10′E海域。随着海温的逐渐降低，1月和2月黄海南部小黄鱼资源继续向外海迁移，即集中在32°20′～34°00′N、123°24′～125°00′E区域，呈现东南-西北向分布，与冬季黄海暖流在黄海南部的主轴位置和流向一致［26］。冬季黄海南部海域受黄海暖流的影响明显，黄海暖流是对马暖流和陆架水混合而成，它主要是在济州岛西侧海域锋区中衍生出来的具有高温高盐特征的海流［26］。因此该海域为黄海南部小黄鱼群体越冬洄游提供了合适的场所。本研究为探寻黄海南部小黄鱼越冬群体提供了重要的空间信息，据此可以通过保护一部分小黄鱼成鱼以维持小黄鱼的合理种群补充水平。

海洋水体环境是海洋鱼类赖以生存的空间，因此可通过对海洋水体环境要素的信息获取、越冬小黄鱼资源的时空分布分析研究小黄鱼越冬的环境偏好，进而可为小黄鱼越冬渔场海况预报提供重要支撑。GAM建模结果显示，黄海南部小黄鱼越冬洄游主要受水深和底层温度影响。随着水深的增加，小黄鱼资源丰度呈现先减少后增加的趋势（图5）。其中水深在30～40 m范围内，小黄鱼资源丰度随水深增加而减少；在40～70 m水深范围内，小黄鱼相对资源丰度总体呈现随着水深增加而明显增加趋势。这与陈新军［27］研究结果基本一致，即越冬期的小黄鱼在黄海区域主要分布在55～75 m。海水盐度的显著变化是支配鱼类行为的一个重要因素，对鱼类渗透压甚至整个鱼类生活史产生重要影响［28］。盐度与鱼类行动的关系主要表现在间接方面，通过水团、海流等起到了影响作用，如暖水性鱼类随着暖流的高温高盐进行洄游［27］。本研究拟合结果显示，底层盐度范围在32.8～33.5时，小黄鱼资源丰度随着盐度的增加而明显增加；当盐度大于33.5时，小黄鱼资源丰度一直较为稳定（图5）。海水盐度变化不是孤立的物理变化现象，它会随着水系、水流等运动而变化，这对于探测小黄鱼越冬分布具有一定意义，后期可通过了解盐度分布及其与海流运动的关系来推测小黄鱼渔场的可能位置。海洋生物受海水温度变化刺激所产生的行为是主动选择最适的温度环境，以使其体温维持在一定的温度范围内。因此海水温度直接影响和控制着海洋生物种群分布及其洄游和繁殖过程，是影响鱼类生态习性的主要环境因子之一［29］。小黄鱼属于暖温性洄游鱼类，其生活习性受温度等水文环境的影响显著［19］，越冬洄游主要是由于环境温度降低所引起的，因此鱼群洄游到渔场时间、资源丰度以及渔期长短均与渔场水温有密切关系。近10年数据拟合结果表明，底层温度对小黄鱼越冬资源丰度影响波动情况较大，但整体上资源丰度随着温度增加而振荡增加（图5），这是因为属暖温性的小黄鱼越冬洄游倾向于水温较高的海域。海洋遥感反演叶绿素a浓度在一定程度上代表该海域的初级生产力水平的高低，常用于估算海洋浮游生物等鱼类饵料的时空分布，因此叶绿素a浓度对小黄鱼资源分布也会产生重要影响。研究结果显示，叶绿素a浓度范围在1.0～2.3 mg·m-3时，小黄鱼资源丰度随着叶绿素a浓度增加而明显增加；当叶绿素a浓度大于2.3 mg·m-3时，小黄鱼资源丰度随着叶绿素a浓度增加而缓慢减少，表明黄海南部小黄鱼越冬洄游资源丰度与叶绿素a浓度并非是完全正相关关系。这可能与在越冬期间鱼类通常减少或停止摄食、而主要依靠索饵期间体内所积累的营养来供应机体能量消耗有关［27］。所以在越冬洄游时期，叶绿素a及饵料的分布和变动在一定程度上难以支配小黄鱼越冬洄游行为。本文主要利用叶绿素a浓度数据进行GAM建模，叶绿素a浓度与其他环境因子相比，其解释率仅为4.6%，即对黄海南部小黄鱼越冬洄游产生影响较小。但后期可通过对比分析索饵期与越冬期海水叶绿素浓度时空动态特征以及此背景下小黄鱼资源丰度情况，深入探究海水叶绿素浓度与小黄鱼索饵和越冬洄游资源丰度之间的时空关系，为预测小黄鱼秋季索饵渔场、越冬渔场位置和制定相关渔业管理政策提供重要参考。

本文建立GAM模型时加入时间（年份）因子作为初始解释变量，可为探究黄海南部越冬小黄鱼资源丰度的年际变动情况、分析小黄鱼越冬适宜生境动态变化提供科学依据。气候环境变化和人类捕鱼活动差异可造成鱼类资源的时空分布变化［30］，并改变其生活史特征和迁移模式［31-32］。探究海洋环境变化与鱼类行为的关系对深入了解和掌握海洋鱼类生活史过程、开发和利用鱼类资源极为重要。所以后期可基于小黄鱼洄游群体资源丰度与环境因子关系，结合海洋卫星环境监测数据对小黄鱼越冬洄游海域和时间、渔期长短等进行预测。