北部湾二长棘犁齿鲷时空分布及其与环境因子的关系

侯 刚，冯钰婷，陈妍颖，王锦润，王思进，赵 辉

（1.广东海洋大学水产学院，广东 湛江 524088；2.广东海洋大学化学与环境学院，广东 湛江 524088）

二长棘犁齿鲷（Evynnis cardinalis Lacèpede,1802）曾名二长棘鲷（Parargyrops edita Tanaka,1916），隶属于鲷科（Sparidae）犁齿鲷属（Evynnis），是北部湾重要经济鱼类。20世纪60―70年代，曾居底拖网渔业产量首位，占比最高达29.3%；但随着过度捕捞加剧和渔场环境改变，北部湾二长棘犁齿鲷资源量严重下降，到20世纪80―90年代几不能形成渔汛[1]。1999年实行休渔制度后，二长棘犁齿鲷曾一度出现“旺发”景象，但近年来该鱼种渔业资源又进入衰退状态。

目前，对北部湾海域二长棘犁齿鲷的研究主要为资源分布与变动[1-3]、种群参数估计[4-7]、种群动力学特征[8-9]、摄食习性[10]、分子标记[11]、开发保护策略[12-13]和休渔效果[14]等，但其资源分布与环境因子关系的研究未见报道。资源分布与海洋环境因子的关系是研究北部湾二长棘犁齿鲷资源变动及资源评估的基础。广义可加模型（generalized additive model，GAM）可很好地拟合渔业资源单位捕捞努力量渔获量（CPUE）与影响因子的非线性关系，已应用于秋刀鱼（Cololabis saira,Brevoort 1856）[15]、拟锥齿鲨(Pseudocarcharias kamoharai,Matsubara 1936)[16]、日本鲭(Scomber japonicus,Houttuyn 1782)[17]和南极鳞虾(Euphausia superba,Dana,1892)[18]等远洋渔业种类，以及小黄鱼(Larimichthys polyactis,Bleeker 1877)[19]等近海种类。本研究根据2013―2014年北部湾渔业资源调查数据，结合遥感海洋环境因子数据，利用GAM模型分析北部湾二长棘犁齿鲷资源与环境因子的关系，为北部湾二长棘犁齿鲷渔业资源开发管理提供参考。

1 材料与方法

1.1 数据来源

为2013―2014年北部湾渔业资源调查数据。调查船为“北渔60011”，船长36.8 m，宽6.8 m，吨位350 t，主机功率661 kW，底拖网上纲长38.5 m，网囊网目 (内径) 2 cm。共4个航次，分别为2013年11月（秋季）、2014年2月（冬季）、2014年5月（春季）和2014年8月（夏季）。以北部湾中东部为调查范围，共设25个站位（图1）。平均拖速3 kn，每站拖曳2 h。渔获全部鉴定至种，逐一计数和称重。资源密度（kg/km2）参照陈作志等研究[1]，按照扫海面积法估算。

图1 北部湾二长棘犁齿鲷调查站位Fig.1 Sampling stations of E.cardinalis in Beibu Gulf

1.2 因子筛选

选取二长棘犁齿鲷的资源密度作为响应变量，海表水温 (Surface sea water temperature,SST)、海表盐度（Sea surface salinity,SSS）、海表叶绿素a浓度 (Sea surface chlorophyll a concentration,Chl-a)、海底地形的数字高程模型（Digital Elevation Model,DEM）、海表面异常（Sea Level Anomaly,SLA）、月份 (Month)、纬度 (Lat)、经度 (Lon) 和离岸距离（Distance）等作为初始环境变量。其中，SST和Chl-a 数据来自美国NASA 的MODIS Aqua 卫星数据 (https://ocancolor.gsfc.nasa.gov/)，数据的时间分辨率为8 d，空间分辨率为4 km。SSS 数据来自Copernicus Marine Environment Management Service(CMEMS) 的 Global Ocean Physical Reanalysis Product 数据（http://marine.copernicus.eu），时间分辨率为月，空间分辨率为1/12°×1/12°。DEM 数据源自Google Earth 高程数据，高程分18 级，空间分辨率为8.85 m。SLA 来自Archiving,Validation,and Interpretation of Satellite Oceanographic (AVISO)数据 (https://www.aviso.altimetry.fr)，时间分辨率为日，空间分辨率为0.25°×0.25°。应用MATLAB 软件提取研究区域SST、Chl-a、SSS 和SLA 卫星遥感数据，去除空值，利用Arcgis10.3 软件将不同空间分辨率的环境数据（SST、Chl-a浓度、SSS和SLA）做重采样，月平均换算成分辨率为0.5°×0.5°的环境数据。利用方差膨胀因子（VIF）度量多重共线性，并对环境变量进行多重共线性检验，筛选适合加入模型的因子。VIF 阈值≥10 时，表明变量间存在严重的多重共线性，则该自变量在建模前予以去除[20]。

1.3 GAM模型构建分析

利用GAM分析二长棘犁齿鲷资源密度和选取因子。GAM的一般表达式[21]：

式中，Y是二长棘犁齿鲷资源密度（kg/km2）；xj为解释变量，即各站位的时空和环境因子；α是适合函数的截距，ε 为残差，ƒi(xj)表示自变量的任意单变量函数，为样条平滑函数。

本研究构建的GAM基础模型：

式中，Y 为二长棘犁齿鲷资源密度，由于二长棘犁齿鲷资源密度存在少量0值，故先Y+1后再进行对数变换；s为自然样条平滑函数；t 为海表温度SST；φ为纬度；A为海表面异常SLA；ρ(Chl-a)为叶绿素a质量浓度；S为海表盐度SSS；λ为经度；d为水深；l 为离岸距离；ε为模型误差，其服从高斯分布。应用R软件中mgcv包实现模型构建和检验[22-23]。利用逐步回归法（Stepwise Method）通过赤池信息准则（Akaike Information Criterion，AIC）和显著性（P值）筛选对模型影响显著的变量，通过F检验及卡方检验评估预测变量加入对模型解释程度的影响[24-25]，确定GAM表达形式。

1.4 模型效果评价

利用AIC 和偏差解释率进行模型选择检验[25]。AIC 值越小，偏差解释率越高，模型效果越好。基于逐步回归法，在AIC 值最小的因子预测函数上依次加入其他因子，得到AIC 值最小的多因子预测模型。当AIC 值不能再减小，构建过程结束，即得到最适模型。AIC 的计算公式如下[25]：

式中，k 为参数数量；L 是似然函数。以上模型构建过程均在R 3.5.2 软件 (R core team) 中实现，其中GAM 由mgcv 包构建[26]。基于Shapiro 等[27]提出的W ′ 正态检验方法，计算统计量，对GAM 模型残差进行正态分布检验。

1.5 渔场重心分析

参考渔场重心分析法，计算各季度二长棘犁齿鲷的资源密度分布重心，分析时空变化特征[28]。渔场重心计算公式：

式中，X、Y 分别为重心位置的经度和纬度；Ci为调查站位i 的资源密度；Xi、Yi分别为调查站位i 的中心经纬度位置；n 为调查站位的总数。

2 结果

2.1 二长棘犁齿鲷的资源分布

北部湾二长棘犁齿鲷分布广泛，出现频率较高，各季调查中该鱼种仅部分站位（秋季站12、冬季站3和站6、春季站19）未渔获。资源密度呈明显的季节性波动，表现为在时间分布上，春季资源密度平均值最高（250.47 kg/km2），最高值为1 578.47 kg/km2；夏季次之（102.10 kg/km2），冬季最低（10.46 kg/km2，图2）。

图2 4 个季度资源密度盒式图Fig.2 Boxplot of stock density by seasons

空间分布上，不同季节二长棘鲷资源密度分布有明显变化，总体为北部高，南部低。秋、冬季资源密度较小，空间分布变化不明显。春、夏季自东北向西南方向减少，春季高密度主要集中在北部湾中部和东北部海域，以20° ― 21° N海域为最高，有4个站位超过500 kg/km2；夏季东北部资源密度较大，但西南部靠近沿岸一带有增加的趋势（图3）。水深梯度分布方面，北部湾二长棘犁齿鲷主要分布在水深11～ 65 m区域，以涠洲岛西南部19～ 50 m的水深梯度急剧变化海域最高，资源密度达702.41～1 578.47 kg/km2（图2、3）。

图3 北部湾二长棘犁齿鲷资源密度时空分布Fig.3 Spatial-temporal distribution of stock density of Evynnis cardinalis in Beibu Gulf

2.2 二长棘犁齿鲷资源密度的GAM模型分析

统计量W ′ 为0.990，达到显著性95%置信水平，即GAM残差符合正态性检验，服从正态分布。通过方差膨胀因子（VIF）对所选8个环境因子进行多重共线性检验，结果见表1。由表1可见，初步设定的环境因子中均未出现共线性问题，故表1中的环境因子均可作为初始变量进行分析。经筛选、验证，该模型选取的时空和环境因子变量包括海表温度（SST，t）、纬度（Lat，φ）、海表面异常（SLA，A）和叶绿素（Chl-a）浓度，得到GAM的表达式为

表1 基于方差膨胀因子的空间和环境因子多重共线性检验Table 1 Multicollinearity test of spatial and environmental factors based on variance inflation factors

利用GAM 模型拟合时空和环境因子变量对二长棘犁齿鲷密度的累计偏差解释率（表2）。GAM模型对二长棘犁齿鲷密度的累积解释偏差为51.6%，调整后R2为0.460。

表2 GAM 模型中已选择变量模拟结果和偏差分析Table 2 Selected variables and deviance analysis for the general additive models

GAM 模型拟合中，通过逐步回归法选取4 个因子，对二长棘犁齿鲷密度影响最大的因子为SST，贡献率为25.9%；其次为Lat、SLA 和Chl-a 浓度，其贡献率分别为11.3%、9.4%和5.0%（表3）。F 检验表明，SST、Lat 和SLA 与二长棘犁齿鲷密度呈显著相关（P <0.05），SST 的显著性最强；Chl-a浓度相关性不显著（P >0.05）。卡方检验表明，SST的检验值最小，非参数平滑效果最好；其次为Lat和SLA，而Chl-a 浓度的检验值不显著（P >0.05），非参数平滑效果低于其他变量。

表3 GAM 模拟空间和环境变量贡献率Table 3 Contributions of the selected spatial and environmental variables in GAMs

北部湾二长棘犁齿鲷资源密度与各预测变量间的效应关系中，环境因子为SST、Chl-a 浓度和SLA，累积贡献率为40.3%（表3）。GAM 分析中，置信区间较窄表明相关性较高。SST 与资源密度呈直线正相关关系（图4_a）。Chl-a 质量浓度在0～ 3.2 mg/m3的范围内与资源密度呈负相关关系，在3.2～5.3 mg/m3的范围内则呈正相关关系（图4_b）。SLA在– 0.19～– 0.1 m 和0.15～ 0.31 m 范围内与资源密度呈负相关关系，在– 0.1～ 0.15 m 范围内则为正相关关系（图4_d）。空间因子Lat 在18.3―21.3° N的范围内与资源密度呈正相关关系（图4_c）。

图4 时空和环境因子对资源密度影响的GAM 分析Fig.4 Generalized additive models (GAMs) for the effects of the spatiotemporal and environmental factors on the stock density

2.3 二长棘犁齿鲷资源密度与环境因子的关系

北部湾二长棘犁齿鲷资源密度与SST、Chl-a浓度和SLA 的叠加分析如图5– 7 所示。二长棘犁齿鲷各季节高值区域，春季主要分布在SST为27.91～ 28.68 ℃、ρ(Chl-a)为0.175～ 2.913 mg•m-3和SLA为– 0.01～ 0.01 m 的海域；夏季主要分布在SST 30.00～ 30.88 ℃、ρ(Chl-a) 0.250～ 0.509 mg•m-3和SLA– 0.19～– 0.06 m 的海域；秋季主要分布在SST 24.43～ 26.26℃、ρ(Chl-a) 0.165～ 1.402 mg•m-3和SLA 0.17～ 0.20 m 的海域。

图5 资源分布与海表面温度的关系Fig.5 Relationship between SST and stock density

图6 资源分布与叶绿素a浓度的关系Fig.6 Relationship between Chl-a and stock density

图7 资源密度分布与海表面异常的关系Fig.7 Relationship between SLA and stock density

2.4 二长棘犁齿鲷渔场重心变动

渔场重心分析表明，北部湾二长棘犁齿鲷的渔场重心有明显的季节移动规律。秋季二长棘犁齿鲷资源密度重心在北部湾西南部海域（107.95° E，19.58° N），冬季资源密度重心向东北方向移动至北部湾中部海域（108.35° E，19.89° N），春季继续向东北方向移动至北部湾东北部海域（108.71° E，20.39° N），夏季向西南方向移动至北部湾东南部海域（102.29° E，19.89° N，图8）。

图8 资源密度分布重心变动Fig.8 Changes in the central gravity of stock density

3 讨论

3.1 二长棘犁齿鲷资源密度与时空因子之间的关系

本研究表明，北部湾二长棘犁齿鲷资源密度季节波动明显，各季节资源密度由大到小依次为春季、夏季、秋季、冬季。渔场重心分析表明，二长棘犁齿鲷为北部湾西南部向东北部的季节性洄游。这种季节性波动与洄游，可能与该鱼种的生殖生长周期密切相关[2]。北部湾二长棘犁齿鲷当年性成熟，于12月―3月产卵[5]，春季时补充群体已遭遇捕捞，因此，5月时资源密度最大。8月时，幼鱼群体索饵料育肥扩散至湾内，11月时当年生鱼已性腺发育，开始向北部涠洲岛海域的产卵场移动。北部湾二长棘犁齿鲷生态寿命虽为7 a左右[5]，但是在遭遇严重捕捞压力情况下，种群年龄结构低龄化，基本上以当年生幼鱼和1龄鱼为主[5]。渔场重心分析表明，北部湾二长棘犁齿鲷季节性洄游较为稳定，与陈作志等[2]研究基本一致，表现为冬季产卵群体游向北部湾东北部，春季产卵群体分散到湾内各处，当年生幼体依然在东北部浅海区育肥生长，至夏季和秋季，鱼群扩散到北部湾各海域[1,3]。北部湾二长棘犁齿鲷的季节性洄游，亦与北部湾的海洋水文环境密切相关[2]。

3.2 北部湾二长棘犁齿鲷资源密度的海洋环境效应

海洋环境是海洋生物赖以生存的基础，海洋生物的栖息分布和繁殖均与海洋水文环境有着密切关系，其中，温度、盐度和深度是主要的阻限[29]。本研究筛选了8个因子作为预测变量进行效应关系分析，环境因子中的SST、Chl-a浓度和SLA和空间因子Lat对模型产生了贡献率（累计51.6%），其中关系最为相关的为海表水温（SST）。二长棘犁齿鲷为底层鱼类，其栖息分布与海表水温无直接关系，但是海表温度为鱼类栖息分布的环境因子之一，与中上层鱼类有直接关系[15,17,30]，而底层鱼类主要通过鱼类繁殖育幼等来影响鱼类的资源分布。海表温度等产卵场环境适宜，鱼卵仔鱼成活率增加，利于资源补充，产卵亲体进行繁殖洄游，导致其资源分布发生季节性变化。

叶绿素a浓度不仅表征浮游植物生物量，而且可以作为以浮游动物为食的浮游动物及海洋鱼类资源丰度的指标[29]。本研究中，叶绿素a与资源密度分布关系不显著（P >0.05），但是其在GAM模型中为第2效应的环境因子，且产生贡献率5%。北部湾二长棘犁齿鲷为当年生当年性成熟鱼类[5]，幼体基本经3～ 4个月生长期即进入捕捞补充，且该北部湾种群在过度捕捞下长期处于种群年龄结构低龄化状态，本研究的二长棘犁齿鲷亦多为当年生幼鱼。有研究表明，北部湾二长棘犁齿鲷在幼鱼阶段，体长30～ 59 mm时主要摄食桡足类，体长60 mm之后摄食对象以鱼类为主，主要为麦氏犀鳕（Bregmaceros macclellandii,Thompson 1940）、康氏小公鱼(Stolephorus commersonnii,Lacepède 1803)等小型饵料鱼类[10]。因此，本研究GAM模型中，叶绿素a通过摄食关系和食物网传递对资源密度分布产生了效应关系。较多研究亦表明，叶绿素a与鱼卵和幼鱼均有较为明显的效应关系[17,31-34]。

海表面异常（SLA）可反映海洋锋面、水团等重要的中尺度海洋动力学特征，含有海浪、潮汐、中尺度涡等海洋动力信息，在渔场分析中有重要作用[34]。北部湾是一个半封闭海湾，受到沿岸水、混合水和外海水的影响，同时近岸的潮汐作用较强，水团、水系等水文特征复杂，而海洋表层水团的大范围输送会导致平均海面高度出现正值和负值[2]。本研究表明，海表面异常SLA显著影响北部湾二长棘犁齿鲷资源密度分布，SLA在-0.19～ -0.1 m和0.15～ 0.31 m范围内，与资源密度呈现负相关关系，在-0.1～ 0.15 m范围内，则为正相关关系。与以往究比较，黄鳍金枪鱼（Thunnus albacares,Bonnaterre 1788 ）、大眼金枪鱼（Thunnus obesus,Lowe 1839）和鲣鱼（Katsuwonus pelamis,Linnaeus 1758）的CPUE在海平面较高处数值较大，而热带大西洋拟锥齿鲨在-0.19 m处CPUE最大[16,30]。这一类差异的原因可能与鱼类生活史状态、栖息特性以及调查捕捞方式有关。空间因子中，纬度Lat与资源密度分布显著相关。北部湾北部水浅平缓，从北向南水深逐渐增加，同时北部湾临近琼州海峡海域海底坡度陡峭，复杂的海底地形和生态环境为二长棘犁齿鲷的繁衍栖息提供了理想的生存环境，并推动了二长棘犁齿鲷东北-西南向的季节性洄游。

4 结论

1）北部湾二长棘犁齿鲷资源密度春夏季高于秋冬季，以春季最高。

2）海表温度、叶绿素浓度和海表面异常是影响二长棘犁齿鲷资源分布的最主要环境因子；纬度是影响资源分布最重要的的空间因子。

3）二长棘犁齿鲷存在明显的季节性洄游，资源密度分布重心，秋季到冬春季从北部湾中南部海域向东北方向迁移至琼州海峡以北海域，夏季向西南方向折返。