棘头梅童鱼(Collichthy lucidus)是珠江河口生态系统重要的经济物种, 也是中华白海豚(Sousa chinensis)重要的饵料生物, 被列入国家重点保护经济水生动物资源名录^[1-2]。现有研究表明, 虽然棘头梅童鱼一年中多批次产卵, 可以承受较高的开发率, 但珠江河口棘头梅童鱼仍处于过度捕捞的状态, 种群出现性早熟现象^[3-4]。探究珠江河口棘头梅童鱼的时空分布及其影响因素是阐明资源分布规律以及进行资源评估的基础, 研究结果可对资源可持续利用以及珠江河口生物多样性保护具有重要意义。

渔业资源时空分布与环境因素之间的关系十分复杂, 具有非线性和不可累加的特性^[5]。广义加性模型能够较好的处理响应变量和多个解释变量之间的关系, 且该模型在构建过程中有较高的灵活性, 因此被广泛用于渔业资源分布和环境因子关系之间的定性分析^[6-8]。在渔业科学研究中, 最佳GAM模型变量组合筛选方法主要分为逐步回归法和全子集回归法, 利用不同的方法筛选出的模型变量不同, 且每种方法都存在一定的缺点^[9-10]。逐步回归法在模型评价全面性方面有限制, 全子集回归法虽然可以评价所有可能的模型, 但当模型中预测变量较多时, 其运行效率较低^[11]。本研究根据2019—2020年春、秋、冬3个季节的珠江河口底拖网调查资料, 采用逐步回归法和全子集回归法, 筛选出棘头梅童鱼最佳的GAM模型变量组合。该方法保留对响应变量有较高贡献率的因子, 评价更多可能模型的同时, 提高了模型运算效率, 可为GAM模型最佳环境变量方法筛选提供参考。

1 材料与方法 1.1 数据来源

渔场数据来自于2019—2020年在珠江河口进行的底拖网资源调查(图1), 调查时间为冬季(2019年12月、2020年1月)、春季(2020年4月)、秋季(2020年9月)。调查船船长21.68 m, 型宽5.8 m, 型深3.3 m, 主机总功率为143.0 kW, 网具的上纲长度36 m, 囊网网目为3 cm, 网衣全长50 m。调查期间, 每站拖曳时间30 min, 拖速控制在3节左右。采集的渔获全部带回岸上进行种类鉴定并测定体长、体质量以及性腺成熟度等生物学数据。

环境数据包括现场测定的站位水深(Depth)、底层盐度(bottom salinity, BS)、底层温度(bottom temperature, BT)、底层溶解氧(bottom oxygen demand, BDO)、底层酸碱度(BpH)等, 叶绿素a 浓度(chlorophyll-a concentra­tion, CHL)数据来自美国国家海洋和大气管理局(https://coastwatch.pfeg.noaa.gov/ ), 时空因子包括经度(Longtitude, Long)、纬度(Latitude, Lati)、季节(Season)。

1.2 数据处理

采用扫海面积法估算资源密度, 计算公式为:

$ D=C / a q, $

(1)

式中, D为每个站点的资源密度, C为每个站点单位小时拖网的渔获量[kg/(网·h)], a为网具每小时实际扫海面积[km/(网·h)], q为可捕系数, 取0.5^[12-13]。

偏度(skewness, S_k)和峰度(kurtosis, K)是量测数据正态分布特性的两个指标, 变异系数(C_V)可以用来衡量数据变异程度的特征大小^[14]。其计算分别为:

$ S_k=\frac{\sum_{i=1}^n\left(X_i-\bar{X}\right)^3 / n}{S^3}, $

(2)

$ K=\frac{\sum_{i=1}^n\left(X_i-\bar{X}\right)^4 / n}{S^4}, $

(3)

$ C_V=\frac{S}{\bar{X}}, $

(4)

式中, S表示标准差, X为平均值, n为样本数量, i表示第i个数值。

1.3 广义加性模型

将Lon、Lat、Depth、BS、BT、BDO、BpH、CHL、Season等因子作为解释变量, 资源密度作为响应变量构建珠江河口棘头梅童鱼的GAM模型, 其中季节为分类离散变量, 其他环境变量为连续变量。

表达式为:

$ \begin{aligned} \lg (D)=& \mathrm{s}(\text { Long })+\mathrm{s}(\text { Lati })+\mathrm{s}(\text { Depth })+\mathrm{s}(\mathrm{BS})+\mathrm{s}(\text { BT })+\\ & \mathrm{s}(\text { DO })+\mathrm{s}(\mathrm{BpH})+\mathrm{s}(\mathrm{CHL})+\text { Season }+\varepsilon ,\end{aligned} $

(5)

式中, D为棘头梅童鱼的资源密度, 对其取对数以满足GAM建模对响应变量须符合正态性的要求; s为自然样条平滑函数, ε为模型误差, 符合高斯分布。

1.4 最佳模型因子组合的选择

采用方差膨胀因子(VIF)对解释变量进行多重线性检验, 避免共线因子叠加作用于模型, 降低误删重要环境因子的概率^[7]。当VIF阈值≥10时, 表示变量存在严重的多重共线性。当只有1个解释变量的VIF阈值≥10时, 去除该解释变量, 并对其他解释变量进行共线性检验; 当多个解释变量的VIF阈值≥10时, 把VIF阈值≥10的解释变量逐步加入其他解释变量进行再次共线性检验。逐一构建单个通过VIF检验的因子与资源密度的GAM模型, 并与初始资源密度GAM模型(Lg(D)=null)的AIC(赤池信息准则, Akaike Information Criterion)值作比较。经过AIC筛选出的因子作基于调整后的R² (Adjust R², adjR²)全子集回归, 将包含较大相对方差解释率的变量且adjR²最大的子集设为最佳GAM模型变量组合。

2 结果 2.1 珠江口棘头梅童鱼资源密度分析

各个季节棘头梅童鱼资源密度偏度均大于0, 数据分布整体呈右偏。冬季和秋季的资源密度频数成尖顶峰分布(K>3), 春季资源密度频数成平顶峰分布(K<0), 各季节渔场密度差异均较大(C_V>0)。棘头梅童鱼季节性分布和集聚强度差异较为明显(S²/m>1; S为标准差, m为均值), 其集聚强度依次为: 秋季>春季>冬季(图2)。

棘头梅童鱼资源密度高值区主要分布在伶仃洋海域。相对于春季, 冬季资源密度高值区向河口方向移动, 而秋季资源密度高值区向外海方向移动(图3)。调查期间, 春、冬两季在万山群岛站位处未捕捞到棘头梅童鱼, 即该区域棘头梅童鱼资源密度为0, 而秋季有较高的资源密度。

2.2 共线性检验

本研究选取的9个环境因子中, Season、BpH、Lati、BT的VIF值均≥10, Long、Depth、CHL、BDO、BS的VIF值<10。依据本研究中共线性检验规则, 把Season、BpH、Lati、BT逐步加入到VIF值<10的环境因子中再次共线性检验, 而Lat、BT仍未通过共线性检验, 即移除Lat、BT后, 其余的预测环境因子VIF值均<10(表1)。

2.3 最佳模型因子的选择

为降低相关因素对模型因子选择的影响, 采用单因子模型AIC值与初始GAM模型AIC值作对比(表2)。

各因子的方差解释率差异较大, 其中CHL、BpH对棘头梅童鱼资源的方差贡献率最高(表3)。Lon 和Depth 单因子GAM模型AIC值大于初始GAM模型AIC值, 将剩余的Season、CHL、BDO、BS、BpH等因子进行全子集回归, 挑选包含CHL、BpH且最大调整后R²的变量组合作为本研究最佳GAM模型的选择因子(图4)。

珠江河口棘头梅童鱼最佳GAM 模型包含CHL、BS、BpH、Season这4个环境要素, 其具体表达式为:

其中, 截距为0.46, Season的斜率为0.28, AIC值为2.90。最佳模型对资源密度的累计方差解释率为96.8%, 调整后的R²为0.901(P<0.001)。其中, CHL贡献率最高, 为58.2%(表3)。

对筛选出的GAM模型进行检验, 残差分布基本符合正态分布, 棘头梅童鱼最佳GAM模型符合要求(图5)。

2.4 资源分布与环境因子的关系

分析最优GAM模型中各环境因子的响应特点发现, BpH相较于其他两个因子的置信误差区间较大, 当BpH>9.6时, 棘头梅童鱼资源密度开始大幅度下降(图6a)。CHL在3.1~5 mg·m^–3时, 棘头梅童鱼资源密度随着叶绿素a浓度的上升而下降, 但置信水平逐渐增加。CHL在5~6.2 mg·m^–3范围时, 其资源密度随着CHL值的增加而上升且该部分置信水平较高。当CHL大于6.2 mg·m^–3时, 资源密度呈现先降低然后上升的趋势(图6b)。BS与资源密度呈负相关关系, 但资源密度随着BS变化而波动较小(图6c)。

3 讨论 3.1 棘头梅童鱼资源密度的时空分布

珠江河口棘头梅童鱼资源密度分布具有明显的季节变化特征, 秋季资源密度最高, 冬季次之, 春季最低, 该规律与珠江口的其他鱼类的季节变动趋势基本一致^[15-16]。棘头梅童鱼虽一年多次产卵, 但冬季和春季是其产卵高峰期, 该段时期出现大量的棘头梅童鱼幼体, 体质量较低, 秋季是棘头梅童鱼补充群体的育肥期, 该时期棘头梅童鱼生物量明显增加, 因此秋季棘头梅童鱼资源密度明显高于其他季节。

棘头梅童鱼作为咸淡水种类, 其空间分布主要与咸淡水区域的位置变化有关。通过进一步分析调查数据, 发现冬季出现棘头梅童鱼主要集中在伶仃洋水域和珠江支流西江河口区域(图7), 结合珠江流域河流径流量推测导致该现象的主要原因是珠江河口分别在伶仃洋水域和西江河口区域形成了较高生产力的咸淡水区域。珠江干流冲淡水汇入伶仃洋水域, 进入伶仃洋水域平均径流量约为5 668 m³/s, 约占珠江流域出海总量的53.4%^[17], 西江是珠江流域的最大支流, 其径流量占珠江流域总径流量的28.3%^[18], 潮流和径流促进咸淡水混合, 相对于其他区域, 河口的咸淡水混合强度最大^[19-20]。咸淡水混合区域汇集了大量的营养盐, 促进浮游动植物的生长, 初级生产力较高^[21], 适宜鱼类栖息, 在珠江干流河口和西江河口区域调查站点均具有较高的资源密度。西江流域冲淡水和外海高盐度水在季风影响下分成不同分支汇入不同的方向, 削弱了西江河口冲淡水向南扩散的动量, 西江河口区域咸淡水面积相对于珠江干流河口水域的咸淡水区域小^[22]。此外, 棘头梅童鱼属于短距离洄游以及集聚型鱼类, 不做长距离的洄游, 珠江口区域的棘头梅童鱼属同一群体^[23], 进一步验证了此推测。

棘头梅童鱼渔业资源高值区在秋季明显靠近河口外部区域, 推测秋季是丰水期, 河流径流量大, 珠江口冲淡水向外海区域扩张, 扩大了棘头梅童鱼等咸淡水种类的分布范围。冬季属于枯水期, 外海高盐水区域向河口入侵, 盐淡水区域向河口收缩, 棘头梅童鱼向河口方向移动^{[16, 23]}。

3.2 GAM模型因子筛选

GAM模型因子筛选过程较为灵活, 采用不同的方法可以筛选出不同的最优变量组合, 设置不同光滑函数的k(节点, Knots)值以及修改逐步回归中环境因子的添加顺序均会对最佳环境因子的筛选结果产生影响, 通常会根据数据情况、目标鱼类的生活史资料等因素采取不同的方法筛选最佳模型变量组合^[24]。如在南印度洋长鳍金枪鱼渔获率与水深、温度关系的研究中, 采用了逐步回归法且以AIC值和显著性水平P为准则建立长鳍金枪鱼渔获率的GAM模型^[25]; 在长江口鱼类群落时空分布特征的研究中, 采用全子集回归的形式并以AIC作为筛选标准确定最佳GAM变量组合^[10], 在沙丁鱼等中上层鱼类的研究中, 将GAM 和GLM相结合探究其渔获量与环境之间的关系^[26]。本研究将逐步回归法和全子集回归法相结合, 筛选出季节、叶绿素、底层盐度、底层酸碱度作为构建最佳棘头梅童鱼GAM模型的最佳变量组合, 该方法不仅对资源密度有较高的方差解释率, 而且在评价更多可能模型的同时, 提高了模型的运算效率。无论采取何种方法准则进行环境因子筛选, 都需要与研究背景结合, 以达到统计学“最佳”和生物生态学“最佳”的结合。

3.3 渔场环境要素对棘头梅童鱼资源的影响

海洋生物的栖息地偏好和生活史特征均与海洋水文环境有着密切的关系, 渔业资源的时空分布受盐度、叶绿素a浓度等环境因子影响较大^{[10, 27-28]}。从鱼类行为学分析, 盐度变化是支配鱼类活动的重要因素^[29], 当盐度发生明显变化时, 会改变鱼类体内外的渗透压平衡, 从而影响其生长和摄食活动^[30], 水团和水流会直接影响鱼类的生长速率和生活史各个阶段的开始时间, 进而影响鱼类种群发育^[30-31]。叶绿素a浓度作为浮游植物生物量的主要指标通过摄食关系和食物网传递与渔业资源密度产生效应关系^[7], 其主要影响棘头梅童鱼饵料生物的丰度分布, 影响棘头梅童鱼的食性和分布^[32]。

本研究中底层盐度对棘头梅童鱼资源密度的贡献度小于叶绿素, 且棘头梅童鱼资源密度随着底层盐度的变化波动较小, 其主要原因是棘头梅童鱼盐度适应范围广且珠江口咸淡水区域底层盐度变化幅度较小, 盐度不是影响珠江口棘头梅童鱼资源密度的首要因素。鱼类对水质的酸碱度均有一定的耐受阈值, 过高或过低的酸碱度对生物生长发育、繁殖都会带来不利的影响, 特别是对底层生物影响大于中上层生物^[33-34], 当底层pH大于9.6时, 棘头梅童鱼资源密度快速下降, 而棘头梅童鱼食性以底层或底栖生物为主^[35], 推测pH对棘头梅童鱼时空分布不仅有直接影响, 还会通过影响其饵料分布产生间接影响。

受中国伏季休渔政策的限制, 夏季采样限制, 因而本研究仅采用了春、秋、冬3个季节的调查数据, 虽然满足建模数据的基本要求, 但其时间序列不够长, 存在一定的偶然性和不确定性。除了物理环境因素以外, 非生物因素如饵料丰度也是影响海洋生物的分布的重要驱动因素之一。在今后的研究中, 应增加调查的频次和调查的覆盖范围, 在建模过程中添加非生物因素, 考虑环境变量间的交互作用, 以期得到更全面的结果。