海洋环境要素是控制鱼类种群分布的重要因子，影响鱼类的生存、生长及繁殖等，环境变动影响鱼类的空间分布和集群(Brander, 2010)。基于生态系统的管理是实现渔业资源可持续发展目标的有效方法(Link, 2002)，确认和理解栖息地环境因子如何影响不同空间和时间尺度上鱼类种群的分布模式是其必要组成部分(Martin et al, 2009)。由此，鱼类种群分布与环境关系的研究成为鱼类研究的重要内容之一。

有关鱼类空间分布与环境的关系，研究方法众多，常用的方法有GAM(Generalized additive model)、GLM(Generalized linear model)、QR(Quantile regression)和分类树等。由于GAM模型能解释鱼类分布和环境因子之间的非线性关系，其近年作为生态学空间格局分析的重要方法之一，被广泛用于探究鱼类种群分布的时空特征及其与环境因子关系的研究(陈新军等, 2007; Furey et al, 2013; Grüss et al, 2014; 王浩展等, 2017; 杨胜龙等, 2017)。随着海洋遥感技术的发展，水色、水温遥感数据已被大量用于海洋渔业研究(樊伟, 2004; 丁琪等, 2015; 崔雪森等, 2018)。

鳀(Engraulis japonicus)是一种洄游性小型中上层鱼类，曾是黄、东海生态系统动力学研究的资源关键种(Tang, 2000; Zhao et al, 2003)，在海洋生态系统中起着承上启下的重要作用。研究表明，越冬鳀空间分布与水温关系明显(马绍赛, 1987; Kim et al, 2005; 李峣等, 2007; 牛明香等, 2012; Niu et al, 2017)。近年来，在环境变动和过度捕捞的影响下，鳀资源量发生明显变化。本研究基于卫星遥感资料和海上调查数据，在定性分析的基础上，利用GAM模型，定量研究黄海越冬鳀资源密度分布与环境因子的关系，旨在为环境变动下的渔业管理提供科学依据。

调查数据来自“北斗”号调查船2000~2015年(2000、2002~2006、2008、2010~2011、2015)1月在黄海中南部(121°~125.5°E、32°~37°N)的底拖网调查，“北斗”号渔业调查船的网具参数：网具网目836目×20 cm，网口周长为167.2 m，网具总长度为83.2 m，囊网网目2.4 cm。站位布设以经、纬度0.5°×0.5°为一个渔区，各年间站位设置一致。调查按预设站位昼夜进行，各站计划拖网1 h左右，平均拖速3n mile/h，实际拖网时间根据调查时的海况确定，详细记录每个站位的渔获组成、位置等相关信息。本研究提取经纬度、水深、下网时间和鳀渔获重量数据。资源密度以单位捕捞努力量渔获量(CPUE, catch per unit effort)表征，即每网每小时的捕获重量[g/(h·net)]。

调查对应时段的水温(SST, sea surface temperature)、叶绿素a(Chl-a)浓度数据来自美国航空航天局(NASA)反演的全球卫星遥感三级数据产品，该数据为月平均数据，空间分辨率为9 km (http://oceancolor.gsfc.nasa.gov)。利用IDL读取原始数据提取研究区(118°~128°E，30°~39°N)的水温、叶绿素a浓度，取1个渔区(0.5°×0.5°)范围内的水温、叶绿素a浓度平均值作为该渔区的值。利用SST计算获得相应渔区的海表温度梯度(TGR, temperature gradient)，计算公式如下：

式中，T_i、T_j为渔区i、j的中心点海表温度。

利用ArcGIS地统计分析模块，在进行数据探索的基础上，Kriging插值生成SST、Chl-a和TGR时间系列分布图，并转换成Raster格式；和相应时间段的鳀资源密度分布图叠加，分析越冬鳀资源密度空间分布与环境因子的关系。

广义可加模型(GAM)是广义线性模型(GLM)的扩展，是一种非参数化或半参数化的回归分析方法，能很好地模拟响应变量和一个或多个预测变量之间的关系(Hastie et al, 1990)。本研究采用GAM模型分析各时空和环境因子对鳀资源密度的非线性影响。由于调查中“零”渔获站位的出现，一般对CPUE加上一个常数，本研究中对所有CPUE加上常数1后取对数，构建的GAM模型如下：

ln(CPUE+1)=s(Year)+s(Hour)+s(Longitude)+s(Latitude)+s(SST)+s(Depth)+s(TGR)+s(Chl-a), family=gaussian

式中，Year表示年份，Hour表示下网时间，Longitude表示经度，Latitude表示纬度，SST表示海表温度，Depth表示水深，TGR表示海水温度梯度，Chl-a表示叶绿素a浓度。模型变量选择基于模型因子的显著性(P < 0.05)，模型选择采用赤池信息量准则(AIC, Akaike information criterion)，将解释变量随机逐个加入GAM模型，得到包含不同因子的模型，AIC值最小者为最佳模型(Damalas et al, 2007)。

受其他因素的影响，有些对鱼类分布影响较大的环境变量有可能对CPUE影响不显著，对该类型变量根据实际情况进行处理。

数据分析利用ArcGIS 10.2、SPSS 17.0和S-plus 8.0软件完成。

由图 1A所知，越冬鳀分布的温度范围较广，从2~17℃都有分布，但高资源密度区主要分布在7~13℃的海域内，且集中在暖水舌及以东的区域。图 1B所示，越冬鳀分布在低叶绿素a浓度海域内，主要范围为0.5~2.5 mg/m³。2000~2015年，越冬鳀分布的温度梯度范围略有变化，2005年开始，鳀分布区域内有高温度梯度水团出现，但纵观所有年份，越冬鳀主要分布在低温度梯度海域，分布范围为0.5~4℃(图 1C)。

图 1(Fig. 1)

图 1 黄海越冬鳀资源密度和环境因子SST、Chl-a和TGR的空间分布 Fig.1 Spatial distribution of stock density of wintering anchovy and environmental factors such as SST, Chl-a and TGR in the Yellow Sea

GAM模型的残差分布及其正态Q-Q图(Quantile-Quantile)常用来检验模型分布假设。图 2A显示，模型残差服从正态分布；图 2B的正态Q-Q曲线分布表明，残差值集中于对角线上并呈线性分布，表明残差趋向于服从均值为零的正态分布；二者均表明GAM模型的分布假设比较合适。

图 1(Fig. 1)

图 2 GAM模型残差的频次分布及正态Q-Q图 Fig.2 Distribution and normal Q-Q plot of residuals in GAM model

依次将各个解释变量逐一加入模型中，每一行表明了加入该因子后解释变量对模型的贡献信息。由表 1可知，年份、下网时间、经度、纬度、温度和水深均为显著性变量，对越冬鳀资源密度分布影响显著，卡方检验P < 0.05；而温度梯度虽然非参数化影响不显著[Pr(Chi) > 0.05]，但加入模型后，F检验P < 0.05，且对模型总偏差解释率大于SST，因此，将温度梯度同样作为显著性因子；叶绿素a浓度对越冬鳀资源密度分布影响不显著，卡方检验、F检验均P > 0.05，但模型加入该因子后，AIC值持续变小，且叶绿素a浓度对越冬鳀的分布也有一定的影响，因此，该变量仍保留在模型中。该模型对CPUE的总偏差解释为44.76%，其中，经度和纬度对CPUE的影响最为显著，分别解释了22.76%和12.64%的偏差；年份、温度梯度、下网时间、水深、温度对CPUE的影响次之，叶绿素a浓度对CPUE的影响最小，解释了0.51%的偏差。

表 1(Tab.1)

表 1 GAM模型拟合结果的偏差分析及显著性因子选择 Tab.1 Summary analysis of deviance for generalized additive models (GAM) and the obvious factors selection 模型因子 Factors added 残偏差 Residual deviance 累积解释偏差 Cumulative of deviance explained% AIC检验 AIC F检验 Pr(F) 卡方检验 Pr(Chi) Null 5641 +年份+ Year 5490 2.68 2277 0.001 0.048 +下网时间+ Hour 5393 4.40 2272 0.018 0.049 +经度+ Longitude 4109 27.16 2159 0.000 0.000 +纬度+ Latitude 3396 39.80 2081 0.000 0.000 +温度+ SST 3333 40.92 2075 0.102 0.001 +水深+ Depth 3250 42.39 2067 0.037 0.017 +温度梯度+ TGR 3145 44.25 2055 0.012 0.300 +叶绿素a浓度+ Chl-a 3116 44.76 2053 0.463 0.323 注：Pr(F)指模型本行与上一行方差分析的F检验值，Pr(Chi)表示模型因子非参数化影响的非线性贡献值 Note: Pr(F) is the P-value from an ANOVA F-ratio test between that row and the previous row of the model; Pr(Chi) indicates the score of the test to assess the non-linear contribution of non-parametric effects

表 1 GAM模型拟合结果的偏差分析及显著性因子选择 Tab.1 Summary analysis of deviance for generalized additive models (GAM) and the obvious factors selection

图 3为各因子对越冬鳀资源密度分布的影响。时间效应方面，越冬鳀资源密度年间变化较大(图 3A)，2000~2015年，资源密度总体呈下降趋势；2000~2002年，资源密度迅速下降；2002~2011年，资源密度年间有差异，但变化幅度较小；2011~2015年，资源密度下降明显。由下网时间与越冬鳀资源密度的关系可以看出(图 3B)，0:00~ 07:00资源密度呈下降趋势，之后处在稳定和缓慢上升的阶段，大约19:00左右，资源密度开始明显上升。

图 1(Fig. 1)

图 3 时空和环境因子对黄海越冬鳀资源密度分布的影响 Fig.3 Impacts of spatio-temporal and oceanographic variables on stock density of wintering anchovy derived from the GAM analysis in the Yellow Sea

空间效应方面，资源密度随经度和纬度变化剧烈。越冬鳀主要分布在122°~124.5°E之间，在124°E以西海域，资源密度和经度呈明显的正相关关系；在124°E以东海域，资源密度和经度呈负相关(图 3C)。越冬鳀集中分布在33°~36°N之间，在34.5°N以南的海域，随着纬度升高，资源密度迅速上升；而在34.5°N以北海域，随着纬度升高，资源密度快速下降(图 3D)。

环境效应方面，越冬鳀资源密度主要集中于60~80 m水深范围内，但水深的影响以大约50 m为分界点，在50 m以浅的水域，随着水深的增加，资源密度升高；而在50 m以深的水域，变化趋势则相反(图 3E)。在7~11℃内，越冬鳀集中分布，水温低于9.5℃时，资源密度与其呈正相关；水温高于9.5℃时，资源密度与其呈负相关(图 3F)。越冬鳀主要分布在温度梯度0.7~2.5℃内，资源密度处于平稳阶段，随着温度梯度增大，达到2℃后，资源密度明显下降(图 3G)。在不同范围内，叶绿素a浓度对资源密度的影响不同，二者之间的关系比较复杂，正相关和负相关交替变化，但总体来看，越冬鳀集中分布在叶绿素a浓度1.2~2.3 mg/m³之间(图 3h)。

2000~2015年，越冬鳀资源密度呈下降趋势，其中，2002~2011年资源密度年间波动不大，但2000~ 2002年和2011~2015年均明显下降。20世纪90年代以前，黄海鳀资源尚未进行规模化开发，之后鳀渔业经过快速发展，在连续几年超过百万吨的捕捞压力下，2003年1月调查发现，资源密度出现历史最低记录，资源接近崩溃边缘(赵宪勇, 2006)。本研究GAM模型结果表明，自2000年开始，2003年资源密度处在最低水平，与上述研究结果一致。鱼类资源密度的影响因素众多，但捕捞压力是最主要的因素，决定鳀的资源丰度。黄渤海鳀年捕捞量(图 4)显示，1989~ 1997年，鳀捕捞产量快速上升，之后两年略有下降，但自1999年开始，捕捞量又开始增加，2003年达到最高值；GAM模型结果显示，2000~2003年资源密度逐年下降(图 3)。自2003年开始，捕捞量快速下降，2009年达到最低，2010年基本持平，之后又明显上升(图 4)；2004~2011年期间，鳀资源密度存在年间波段，但总体看来，2010年和2011年呈上升趋势，而2015年却快速下降(图 3)。资源密度与捕捞量的关系充分印证了捕捞压力对鳀资源密度的影响。

图 1(Fig. 1)

图 4 黄渤海鳀年捕捞量(数据来源：中国渔业统计年鉴) Fig.4 Fishing yield per year of anchovy in the Yellow Sea and Bohai Sea (Data source: China Fishery Statistics Yearbook)

GAM结果表明，调查捕获的鳀资源密度一天之内变化明显，白天所捕获的资源密度较夜间稳定，且从白天到夜晚，资源密度不断升高，而从夜晚到白天，资源密度逐渐下降。这是因为冬季受强烈北风的影响，水层垂直混合作用明显，黄海水域基本无温跃层出现(邹娥梅, 2001; 葛人峰等, 2006)。鳀垂直分布昼夜变化明显，在无温跃层水域，白天常分布于中下层或近底层，且易集成小群；夜间则主要分散分布于中层至中上水层(马绍赛, 1989; 赵宪勇, 2006)。林德芳(1997)研究显示，对于鳀的不同集群模式，其集群密度不同，中层型高于近底层型，近底层型高于表层型。因此，白天到夜间，鳀的集群模式由近底层型转变为中层型，资源密度不断升高；夜间到白天，则反之。

黄海暖流对黄海越冬鳀的分布有重要影响，越冬鳀主要密集于黄海暖流入侵所形成的暖水舌边缘(陈聚法, 1997; 李峣等, 2007)，图 1A中鳀的分布与其相符。GAM分析显示，经纬度对越冬鳀资源密度影响最大，贡献率占模型总贡献率的79%。在124°E以西海域，越冬鳀资源密度随经度升高而快速上升；而在124°E以东海域，由于调查数据减少，数据相关性变差。在纬度分布上则以34.5°N为分界线，34.5°N以南的海域，越往北，资源密度越高；而34.5°N以北海域，越往北，资源密度越低，且变化明显。赵宪勇(2006)的研究表明，越冬期间不同大小个体鳀的分布有明显的区域特征，大个体鳀主要分布于123°E以东海域，而小个体则主要分布于调查水域的北部和西部，以西部浅水区为主。本研究结果与其一致。

越冬鳀资源密度分布与水深关系密切，模型贡献率为1.47%。与其他因子不同的是，虽然资源密度主要集中在水深60~80 m范围内，但水深对资源密度的影响却以大约50 m为界。由图 5可知，从浅水区到深水区，资源密度逐渐增大，因此，在浅于50 m的水域内，水深与资源密度呈正相关。随着水深的继续增加，每个调查站位的资源密度没有降低；但由于GAM结果表现的是一种相对效应，随着水深增加，该深度范围内的站位数量减少，因此，总体看来，深于50 m的水域，资源密度随水深增加而降低。

图 1(Fig. 1)

图 5 越冬鳀资源密度与水深等值线分布[g/(h·net)] Fig.5 Map of wintering anchovy stock density and water depth contour [g/(h·net)]

水温对越冬鳀资源密度分布影响显著，越冬鳀密集分布在7~11℃范围内；小于9.5℃时，资源密度随温度升高而增加；大于9.5℃后，资源密度随温度升高而降低。这是因为水温是鱼类生活环境中最基本最重要的因素，其变化对鱼类的生活、生长及饵料的生长均产生影响，这直接或间接地影响着鱼类资源量的分布。鳀作为小型中上层鱼类，对水温变化敏感(李向心, 2007)，海水温度的变化影响着越冬鳀的纬向分布(牛明香等, 2012; 陈云龙, 2014; Niu et al, 2017)。马绍赛(1989)和李峣等(2007)的研究表明，越冬鳀的适温范围为8~15℃，密集分布区为11~13℃，本研究结果略有差异。这是因为鱼类不同生活史阶段适温范围不同(Bellido et al, 2001)，种群大小、年龄结构和捕捞状况等均影响越冬鳀的分布(Loots et al, 2010)。本研究所用为2000~2015年数据，时间跨度大，前期研究数据分别为1985~1987年和2002~2003年；随着捕捞压力的变化，黄海越冬鳀的年龄结构变化明显(赵宪勇, 2006; 张俊, 2011)，因此，适温范围有所差异。

越冬鳀的高资源密度区主要分布在温度梯度相对小的区域内(图 1C)。GAM结果表明，越冬鳀主要分布在大约0.7~2.5℃的温度梯度范围内。在适宜的范围内，随着温度梯度增大，越冬鳀资源密度由基本稳定到明显下降。这是因为越冬鳀主要密集分布于黄海暖流暖水舌边缘(李峣等, 2007)，该区域温度梯度不大；而在温度梯度大的浅水区域，则主要分布着小个体鳀鱼(赵宪勇, 2006)。

越冬鳀资源密度分布与环境因子关系密切，资源密度集中分布在一定的海表温度、叶绿素a浓度和海表温度梯度范围内。GAM定量分析表明，年份、下网时间、经度、纬度、水深、温度和温度梯度对越冬鳀资源密度分布影响显著。2000~2015年，越冬鳀资源密度总体呈下降趋势，年间波动明显，个别年份资源密度呈上升趋势；下网时间对捕获的鳀资源密度影响明显，从白天到夜晚，资源密度逐渐升高，而从夜晚到白天，资源密度不断降低；调查范围内，经度对资源密度的影响以124°E为界，纬度对资源密度的影响则以34.5°N为界，两端呈现相反的变化趋势；50 m以浅的水域，资源密度随水深增加而增加，50 m以深的水域，资源密度随水深增加而降低；越冬鳀集中分布在7~11℃范围内，水温明显影响着资源密度的变化趋势；最适宜的海水温度梯度和叶绿素a浓度分别为0.7~2.5℃和1.2~2.3 mg/m³。