渤海口虾蛄时空分布及其主要影响因子的空间异质效应*

胡贝尔, 韩青鹏, 关丽莎, 徐开达

(1.中国水产科学研究院黄海水产研究所 农业农村部海洋渔业可持续发展重点实验室, 山东省渔业资源与生态环境重点实验室,山东 青岛 266071; 2. 大连海洋大学水产与生命学院, 辽宁 大连116023; 3. 中国海洋大学水产学院, 山东 青岛 266003;4. 浙江省海洋水产研究所 浙江省海洋渔业资源可持续利用技术研究重点实验室, 浙江 舟山 316021)

口虾蛄(Oratosquillaoratoria)隶属甲壳纲(Crustacea)口足目(Stomatopoda)虾蛄科(Squillidae)口虾蛄属(Oratosquilla),是中国、日本、韩国及东南亚沿岸小型渔业的重要捕捞对象之一,其具有重要的经济和生态价值[1-3]。在过度捕捞和近海生境退化等多重压力下,全球口虾蛄资源普遍衰退,捕捞群体以1龄甚至0龄补充群体为主[4-5]。2021年,农业农村部将口虾蛄新增为我国海洋伏季休渔期间专项捕捞品种,规定了具体的作业时间、作业海域、作业船数、作业类型、渔具网目尺寸和捕捞限额等,以期在口虾蛄渔业资源保护和利用两者间达成平衡。在此背景下,预测口虾蛄资源密度的时空分布,识别主控环境因子并解析其对资源分布的影响, 将为口虾蛄渔业的科学精细化管理提供生态学基础。

地理加权回归(Geographically weighted regression, GWR)模型是对传统回归模型的扩展,将全局参数转换为随地理位置变化的局部参数[6]。不同于广义线性模型(Generalized linear models, GLM)和广义可加模型(Generalized additive models, GAM),GWR模型中解释变量(如环境因子)对响应变量(如资源密度)的影响是不固定的,将研究区域分割成不同的子区域,并在每个子区域建立局部回归方程,以分析解释变量对响应变量的局部效应。因此,GWR模型通常具有更高的预测准确性,能够通过反映解释变量的局部效应来分析响应变量的时空分布变化[7]。

此前关于口虾蛄时空分布的研究大多采用传统统计方法,研究环境因子对口虾蛄密度分布的全局效应。例如,吴强等[5]和刘泽修等[8]通过Pearson相关分析分别研究了莱州湾与辽东湾口虾蛄的空间分布及其影响因素;李鹏程等[9]利用GLM方法,研究了山东近海口虾蛄的空间分布特征及其与环境因子的关系。尚未有报道分析渤海口虾蛄密度分布影响因子的局部效应。本研究为分析整个渤海空间尺度的口虾蛄时空分布模式,采用GWR模型,识别影响口虾蛄密度分布的主控环境因子,并解析其对渤海口虾蛄密度分布的局部效应,以期为口虾蛄渔业资源管理提供科学基础。

1 材料与方法

1.1 数据来源

本研究数据来源于中国水产科学研究院黄海水产研究所于2017年春季(4—6月)和秋季(9—11月)在渤海开展的渔业资源底拖网调查。各调查航次的口虾蛄密度分布及底层水温等温线图见图1,调查租用双拖渔船,渔船功率为108 kW。专用调查网具网口高6 m、宽22.6 m、周长为1 740 目, 网目 63 mm, 囊网网目 20 mm。每站拖网1 h,拖速为3 kn。但因天气、生产渔船作业、站位水深和底质限制等原因, 个别站位拖网不足1 h,各月份调查站位也有所不同。每个航次调查站位数具体为:4月36个、5月39个、6月46个、9月48个、10月43个、11月39个。记录所有捕获口虾蛄的数量(ind.)和质量(g)。将各站位质量数据按照1 h进行标准化处理,得到各站位口虾蛄密度(g/h)数据。

图1 2017年春季(4—6月)和秋季(9—11月)渤海底拖网渔业资源调查中口虾蛄密度和底层水温分布图

1.2 统计分析

应用地理加权回归(GWR)模型分析口虾蛄密度分布与环境因子之间的关系,并筛选显著影响因子。其基本表达式为:

logD=log(D+1) ,

(1)

(2)

式中:D为口虾姑密度;xjk为第j个局部区域处第k个自变量的观测值;m为自变量的个数;βj0为局部区域j处的截距参数;βjk为局部区域j第k个自变量的局部回归系数;εj为局部区域j处的随机误差。该模型采用了一个高斯核函数,通过自动带宽选择方法确定了一个最优的自适应带宽(即一个反映局部样本大小的固定量)。

本研究初步考虑了海水底层温度(Bottom temperature, BT)、底层盐度(Bottom salinity, BS)、水深(Depth)及5个底质变量等8个环境因子对口虾蛄出现概率的影响。底质变量分别为含砂量等级(Order in the sand percentage, Ord_Psand)、黏土含量等级(Order in the sand percentage, Ord_Pclay)、底质粒径分布的平均粒径等级(Order in mean grain size, Ord_Mz)、偏度等级(Order in skewness, Ord_Sk)和峰度等级(Order in kurtosis, Ord_kt)。底质数据来源于文献[10]的渤海表层底质类型和粒径分布图(数据见表1和表2)。因该底质粒径分布图为等值线图,故本文所提取的底质数据为整数型,表示各底质因子的大小等级[10]。

表2 底质变量的等级与相应范围

备选底质因子间存在较强的共线性(见表3),但不通过模型选择很难确定应该保留哪几个因子。因此,采用逐步回归法[11],根据赤池信息准则AICc筛选显著影响口虾蛄出现概率的主要环境因子和最优模型。另外,采用交叉验证法分别验证春季和秋季最优模型的预测准确性。交叉验证的具体步骤:(1)随机抽样80%作为训练数据,剩下20%为数据;(2)用最优模型拟合训练数据,并根据实验数据的环境因子预测实验站位的口虾蛄密度;(3)通过线性回归(y=bx)分析实验站位的观测密度与预测密度的关系,并记录回归系数和模型解释率R2;(4)重复上述过程500次,分析500个回归系数和R2的中值和分布。利用最优模型预测春季和秋季各月份渤海口虾蛄出现概率的空间分布。

表3 底质变量间的相关性系数

2 结果

2.1 最优模型和主要影响因子

春季,最优GWR模型(M25)包括四个影响因子,按贡献率从大到小依次为BT、Ord__Kt、Depth和Ord_Mz(见图2(a))。该模型的累积偏差解释率为50.4%,AICc值为538.7,继续增加影响因子会导致AICc值增加(见图2(c))。秋季,最优GWR模型(M21)包括三个影响因子,按贡献率从大到小依次为Ord_Mz、Depth和 Ord_Kt(见图2(b))。该模型的累积偏差解释率相对较低,为28.3%,AICc值为544.6,继续增加影响因子会导致AICc值明显增加(见图2(c))。交叉验证结果表明,春季和秋季最优模型均具有较高的预测准确性,斜率接近于1且R2值超过0.9(见图3)。

(图(c)中模型对应的最小AICc值以红色标注。Minimum AICc of the optimal model was marked in red in panel (c).)

(垂直虚线为R2的均值。The vertical dotted line is the mean of R2.)

2.2 环境因子影响的空间异质性

春季,GWR模型中各环境因子对口虾蛄密度分布的局部效应存在明显的空间异质性(见图4)。底层温度BT的局部效应系数均为正值(0.8～2.0),即口虾蛄密度随着BT升高而增加;其相对高值区(1.4～2.0)主要位于莱州湾和辽东湾(见图4(a)),低值区主要位于渤海湾及邻近水域(见图4(a))。底质粒径分布的峰度等级Ord_Kt的局部效应系数正负不一,其对口虾蛄密度分布的影响在辽东湾和莱州湾表现为正效应(即底质粒径分布的峰度等级越大,口虾蛄密度越高),而在渤海中部呈较强的负影响效应(见图4(b))。水深Depth对口虾蛄密度分布的影响除在辽东湾表现为负效应,在渤海其他水域主要表现为正效应(即口虾蛄密度随水深增加而增加)(见图4(c))。平均粒径等级Ord_Mz对口虾蛄密度分布的正效应主要位于辽东湾及莱州湾湾口东北部,负效应主要位于莱州湾和渤海中部(见图4(d))。

图4 春季最优GWR模型中各环境因子对口虾蛄密度的效应系数的空间分布

在秋季,环境因子对口虾蛄密度分布的影响也存在明显的空间异质性(见图5)。其中,平均粒径等级Ord_Mz对口虾蛄密度分布的影响为正效应,即Ord_Mz越大(肉眼观测的底质粒径越小)口虾蛄密度越高;在莱州湾附近,Ord_Mz的局部效应系数明显更大,影响也更大(见图5(a))。水深对口虾蛄密度分布的影响主要表现为正效应,即口虾蛄密度随水深增加而增加;在莱州湾和渤海湾,水深因子的正效应相对较强;另外,水深对口虾蛄密度分布的影响在辽东湾湾口(水深较深,多在25 m以深)表现为较弱的负效应(图5(b))。底质粒径分布的峰度等级Ord_Kt对秋季口虾蛄密度分布的影响主要表现为正效应,在莱州湾其局部效应明显更强,即口虾蛄密度随Ord_Kt增大而增加;在其他区域,Ord_Kt的局部效应系数基本在-0.5～0.5之间,对口虾蛄密度分布的影响相对较弱(见图5(c))。

图5 秋季最优GWR模型中各环境因子对口虾蛄密度的效应系数的空间分布

2.3 密度分布

基于GWR模型预测,在4和5月渤海口虾蛄高密度分布区主要位于渤海湾和莱州湾及邻近水域,在辽东湾及邻近渤海中部水域口虾蛄密度相对较低;这一预测结果与实际观测值基本相符(见图6(a)和6(b))。2017年6月,渤海口虾蛄广泛分布于整个调查区域,且密度明显高于4月和5月,主要高密度分布区仍位于渤海湾和莱州湾以及邻近水域,另外辽东湾也出现相对密度高值区(见图6(c))。

(图(b)标注2021年5月伏季休渔期间专项特许捕捞口虾蛄作业区域(即多边形内部区域)。The open areas of special licensed fishing for mantis shrimp are marked with polygons in panel (b) during the fishing moratorium in May, 2021.)

秋季,因最优GWR模型中所包含的环境因子均为静态因子,9—11月的预测密度分布一致。秋季,口虾蛄密度高值区集中分布于渤海湾及邻近水域和包含莱州湾西北部在内的渤海中部沿岸水域,在辽东湾存在散布的高密度区;在渤海中部深水区存在明显的条带状低密度分布区,这与该区域口虾蛄较低的密度观测值相一致(见图6(d)—(f))。

3 讨论

3.1 口虾蛄空间分布的季节变化

渤海口虾蛄始终生活在渤海海域,12月至次年3月在渤海中部越冬,3月下旬开始向近岸产卵区洄游[12]。本文结果显示4月和5月口虾蛄高密度分布区主要位于渤海湾和莱州湾及其相近区域,渤海中部密度较低。辽东湾相对较低是由于辽东湾北部近岸区域未设调查站位,董婧等[13]发现,4—5月辽东湾北部近岸区域分布着数量较多的口虾蛄;产卵盛期(5月)后,6月口虾蛄开始在渤海广泛分布,进行索饵。本研究中观察到6月口虾蛄密度明显高于4和5月,可能是因为口虾蛄开始离开近岸产卵场。秋季属于索饵期,本文研究发现口虾蛄主要分布于渤海湾及邻近水域和包含莱州湾西北部在内的渤海中部沿岸水域,9和10月的预测结果与拖网调查中站位密度有较好的契合度。从拖网调查中站位密度可以看出,秋季随着水温的降低,口虾蛄开始向渤海中部聚集,准备越冬。因此,包含静态因子的秋季最优模型可以较好地预测9、10月口虾蛄密度分布,但对11月口虾蛄分布的预测准确性降低。 休渔期后,由于捕捞强度的影响,口虾蛄密度会有所下降[9,12]。因此,本研究拖网调查结果和模型预测结果均显示秋季期间,口虾蛄密度相对下降。受生活史特征的影响,渤海口虾蛄形成春季近岸密度较高、秋季远岸密度较高的规律。

3.2 口虾蛄密度分布的主要影响因子及其空间异质效应

在生活史的不同阶段,物种对环境因子的需求有所差异,各阶段主控环境因子的空间异质性导致了物种分布的差异[14]。温度作为影响海洋生物分布的重要生境因子,可直接影响口虾蛄新陈代谢等生命活动[15],其对口虾蛄分布的影响在春季最显著。春季是口虾蛄的主要繁殖时期[16],4月开始游向近岸产卵,5月为繁殖盛期,口虾蛄聚集在水温12～18 ℃的近岸浅水区。在本文中,模型预测的春季口虾蛄密度随着温度的升高而增加,这与春季渤海最高海底温度均未超过口虾蛄分布的适宜温度阈值上限有关。春季,口虾蛄在莱州湾及邻近水域分布较密集,可能是因为该区域水温相对较高,在口虾蛄的适宜温度范围内;且该区域水温分布梯度明显,更能反映温度变化对口虾蛄密度分布的影响,因此水温的局部效应系数较高。相反,在渤海湾水温对口虾蛄密度分布的局部正效应相对较弱,这主要是因为该区域4—6月的水温接近口虾蛄的适宜温度范围,且水温分布较为均匀(见图1)。由此可知,渤海底层水温分布的空间异质性导致水温对口虾蛄密度分布的局部效应存在明显的空间差异。使用广义线性模型(GLM)和广义加性模型(GAM)等全局模型难以捕捉到这种空间异质性,进而影响了口虾蛄密度分布预测的准确性。

在春、秋两季,水深对口虾蛄密度分布的影响均以正效应为主,负效应出现在辽东湾及其湾口以北等水深较深区域(多在25 m以深)。这意味着,渤海口虾蛄更喜分布于特定的水深范围;在25 m以浅水域,渤海口虾蛄更喜分布于较深水域;当水深超过25 m,口虾蛄的密度随水深增加而减少。另外,Zhang等[17]利用GAM模型对海州湾及邻近水域口虾蛄栖息地适宜性的研究表明,其适宜性随水深增加而降低。此研究结果的差异可能是由不同口虾蛄地理种群的局部适应性造成的,也可能是因为它们分别基于不同的模型技术(即全局模型GAM和地理空间模型GWR)。

在春、秋两季,平均粒径等级Ord_Mz均为影响口虾蛄密度分布的主要影响因子,但其对口虾蛄分布的影响效应存在季节差异。春季,在平均粒径等级相对较低(Ord_Mz为2～4,平均粒径约为3～6 Φ,底质类型主要为极细砂、粗粉砂和中粉砂)的辽东湾和莱州湾湾口(见表1、2和文献[10]的图4)口虾蛄密度随Ord_Mz增大而增加;换言之,口虾蛄更集中分布于平均粒径等级为4、平均粒径接近6 Φ的中粉砂底质类型。相反,在平均粒径等级相对较高(Ord_Mz为4～5,平均粒径约为5～7 Φ,底质类型主要为中粉砂和细粉砂)的莱州湾和渤海中部(见表1、2和文献[10]的图4),Ord_Mz具有明显的负效应;这意味着,口虾蛄更集中分布于平均粒径等级为4、平均粒径为5～6 Φ的中粉砂底质类型。无论在春季还是秋季,Ord_Mz的效应在渤海湾接近于0或正效应较弱,主要是因为此区平均粒径等级比较均匀(见文献[10]的图4)。秋季,在渤海湾中,平均粒径等级Ord_Mz对口虾蛄密度分布的影响为正效应,这意味着口虾蛄更偏爱平均粒径等级高的底质类型。李鹏程等[9]研究发现海州湾口虾蛄主要集中于淤泥质粉砂底质;许莉莉等[18]发现海州湾口虾蛄主要栖息在砂、砂-粉砂-黏土和黏土质砂等类型的底质中,这与本文研究结果较为一致。

底质粒径分布的峰度等级Ord_Kt是影响口虾蛄密度分布的另一主控环境因子。在莱州湾,其对春、秋季口虾蛄密度分布的影响均表现为明显的正效应(见图4(b)和图5(c));结合渤海底质粒径分布的峰度分布图(文献[10]的图4d),发现当Ord_Kt较低(即峰度小于3)时口虾蛄密度随峰度增大而增加。在渤海中部Ord_Kt相对较高(即峰度大于3)水域,Ord_Kt对口虾蛄密度分布具有稳定的负效应(见图5(c)),即口虾蛄密度随峰度的增大而减小。概括而言,口虾蛄偏爱的峰度范围应在3左右。在辽东湾,底质粒径分布的峰度存在较高空间异质性(见文献[10]的图4d),其对口虾蛄密度分布的影响基本表现为正效应,但效应强度在春季较高。粒径分布的平均粒径、偏度和峰度直接影响底质的质地,进而影响其持水能力[19-20]。由于水的比热容明显高于固态底质,底质的持水量直接影响其热容,进而可影响穴居生物的体温调节[20-21]。因此,在研究穴居海洋生物密度分布和栖息地适宜性时,建议结合其生物学和生态学以更为科学合理地筛选主控环境因子。