郑好好,杨晓明,2,3,朱江峰,2,3
1.上海海洋大学 海洋科学学院,上海 201306
2.国家远洋渔业工程技术研究中心/大洋渔业资源可持续开发教育部重点实验室,上海 201306
3.农业农村部大洋渔业开发重点实验室/农业农村部大洋渔业资源环境科学观测实验站,上海 201306
中西太平洋鲣 (Katsuwonuspelamis) 是世界上重要的渔业资源,为全球人口供应动物性蛋白质,并在出口创汇、渔业转型升级及远洋渔业相关产业发展中作出了重要贡献。金枪鱼围网渔业在我国中西太平洋远洋渔业的地位举足轻重,而鲣是围网渔业中的重要捕捞对象[1]。2020 年中西太平洋海域捕捞量的72%为金枪鱼围网捕捞,其中鲣占围网总捕捞量的82%[2]。多数研究认为,在影响鲣资源的重要因素中海洋环境因素不可忽视[3-5],因此,开展有关海洋环境与鲣资源关系的研究具有重要意义。
目前一些研究方法被广泛应用于探讨鲣资源空间分布与海洋环境之间的关系,如随机森林法则[6](Random Forest,RF)、BP 神经网络[7-8],广义加性模型[9](Generalized Additive Model,GAM) 和最大熵模型[10-11](Maximum Entropy Model,MaxEnt)等。然而,这些模型并未考虑海洋环境因子影响的空间异质性问题。针对传统线性回归模型所忽略的空间异质性问题,地理加权回归模型[12-13](Geographically Weighted Regression,GWR) 在一定程度上有所改进,但该模型只能反映各变量的平均尺度,未考虑不同变量的空间异质性尺度差异,无法体现不同环境因子的多尺度效应,因此也存在一定的估计偏差。
近年来,在分析影响因素空间异质性的实证研究中,逐渐纳入多尺度地理加权回归模型[14](Multi-scale Geographically Weighted Regression,MGWR),随着相关的统计推断不断地补充完善,MGWR 模型已较普遍地应用于研究中[15-16],且已证明其在空间尺度差异和异质性研究上有较好的拟合效果[17]。MGWR 模型考虑了各影响因子的尺度差异,本研究采用该模型方法,运用到中西太平洋鲣渔获率与海洋环境因子关系研究中,同时为探究MGWR 模型的精度,选取了GAM 和GWR 模型为对照,对比分析这3 种模型的拟合优度,探讨了不同环境因子对中西太平洋鲣渔获率空间分布的影响,选出最合适的模型,为鲣资源的养护管理和合理开发利用以及我国金枪鱼围网渔船的生产提供参考依据。
1.1.1 渔业数据
渔业数据来源于中西太平洋渔业委员会 (Western and Central Pacific Fisheries Commission,WCPFC) 所公布的中西太平洋海域围网渔业捕捞数据,本文选取2005—2019 年140°E—160°W 和15°S—15°N 海域范围内的捕捞对象为鲣的数据进行研究,时间分辨率为月,空间分辨率为1°×1°,包括年份、月份、经纬度、捕捞天数及渔获量等。
1.1.2 环境数据
选取与鲣活动和栖息相关的环境因子,包括不同水层 (5、55 、100 、150、200 m) 的温度、盐度、东西向和南北向海水流速,以及净初级生产力。所有环境因子均来自参与第六次国际耦合模式比较计划 (Coupled Model Intercomparison Project,CMIP6) (https://esgf-node.llnl.gov/search/cmip6/) 气候-地球系统模式研发的美国国家大气科学研究中心 (National Center for Atmospheric Research,NCAR) 及美国国家大气海洋局地球流体动力学实验室 (National Oceanic and Atmospheric Administration-Geophysical Fluid Dynamics Laboratory,NOAAGFDL),其中2014 年及以前的环境数据用于评估CMIP6 地球系统模式在历史时期的模拟表现,2015 年以后使用了最新的综合评估模型和排放数据[18],本文选取其中的2005—2019 年的环境数据,并通过Matlab R2018b 软件将其与渔业数据进行匹配。时间尺度为月,空间尺度为1°×1°。
1) 计算2005—2019 年1°×1°各单元渔区内的累计渔获量和累计作业时间 (d),从而获得2005—2019 年各单元渔区内渔获率,即名义单位捕捞努力量渔获量 (Catch per unit effort,CPUE)。
2) 考虑到捕捞时的效率容易受各种因素的影响而产生差异,本文采用名义CPUE 表示渔获率的分布,其计算公式为:
式中:i,j表示作业渔船的经纬度位置;YCPUE表示名义单位捕捞努力量渔获量(t·d-1);Ucatch表示对应地点作业渔船的累计渔获量 (t);fdays表示在该作业位置累计作业时间 (d)。
1.3.1 环境因子的选取
当各因子间存在多重共线性时,会出现过拟合而导致模型泛化能力降低。本文采用方差膨胀因子(Variance inflation factor,VIF) 判断各海洋环境变量是否存在多重共线性,提取VIF<7.5 的环境变量[19]。
1.3.2 GAM
GAM 是一种全局回归模型,它使用未指定的(即非参数的)平滑函数来代替线性协变量函数。由于该方法可以有效地分析环境因子和渔获率之间的非线性关系[20-21],因此被广泛运用于渔业研究中。利用GAM 模型建立CPUE 与环境因子之间的关系式为:
式中:S(.)为各环境因子的样条平滑函数;εGAM为模型拟合残差。
1.3.3 GWR 和MGWR
GWR 模型是在传统的全局回归模型的基础上,将数据的地理位置纳入回归参数,在估计局部参数时考虑相邻点的空间权重[22]。其模型表达式如下:
与经典GWR 模型相比,MGWR 模型允许每个变量各自不同的空间平滑水平,这降低了估计的偏误,同时也产生了更真实有用的空间过程模型,MGWR 模型的计算公式如下:
式(3)、式(4)中:xij是变量xj在i点的值;(ui,vi)代表样本点的空间地理位置;βj(ui,vi) 是i点上的第j个回归参数,是地理位置的函数,当j=0 时,β0(ui,vi)为i点的回归常数;k为回归系数的总个数;βbwj代表了第j个影响因子回归系数使用的最佳带宽;εi为随机误差。MGWR 模型的各回归系数βbwj均基于局部回归所得,且带宽具备特异性[23],在经典GWR中,βj所有影响因子的带宽都相同。
MGWR 模型的核函数和带宽选择准则也与经典GWR 模型一致,本文采用Gauss 核函数,带宽选取准则为修正的赤池信息准则 (Corrected Akaike Information Criterion,AICc)。对比可知,MGWR 通过推导出响应变量和不同预测变量之间条件关系的单独带宽,允许不同的过程在不同的空间尺度上运行。MGWR 使用反向拟合算法进行校准,用GWR 参数估计初始化反拟合过程。基于这些初始值,校准过程以迭代的方式工作,在每次迭代中,所有的局部参数估计和最优带宽都被评估。当连续迭代的参数估计的差值收敛于指定阈值时,迭代终止。本研究中,收敛阈值取为10-5。
本文MGWR 模型的计算基于美国亚利桑那州立大学空间分析研究中心 (SPARC) 开发的MGWR 2.2 软件 (https://sgsup.asu.edu/SPARC),该软件在进行模型拟合时可选择对环境变量进行标准化处理,其表达式为:
式中:Zij为标准化后的变量值;Xij为实际变量值;和σ分别为实际数据的均值和标准差。采用Arc-GIS 10.5 软件制作地图。
1.3.4 模型性能评估
为评估模型的性能,本文对比了GAM、GWR和MGWR 模型的AICc、残差平方和 (Residual sum of squares,RSS)、拟合优度 (R²) 和校正后的拟合优度 (AdjustedR²)。较低的AICc 值、RSS 值和较高AdjustedR2值意味着模型具有较好的拟合效果。
通过多重共线性诊断,筛选出以下5 个海洋环境因子,结果见表1。
表1 解释变量间方差膨胀因子Table 1 Variance inflation factor among explanatory variables
GAM、GWR 和MGWR 模型的计算结果参数见表2。与GAM 模型相比,GWR 模型AICc 值略高,调整后的R²显著提升;两个模型的残差平方和 (Residual Sum of Squares,RSS) 对比显示,GWR模型的RSS 显著降低,表明GWR 模型的拟合度高于GAM 模型。在比较GWR 模型和MGWR 模型的参数时,发现MGWR 的RSS 和AICc 值明显较低,显示出MGWR 模型的结果优于GWR 模型。此外,MGWR 模型调整后的R²从0.846 增至0.871,进一步证明了MGWR 模型的拟合优度最好。总体而言,MGWR 模型的性能明显更好,可以用来解释CPUE 的空间分布。
表2 GAM、GWR 和 MGWR 不同回归模型性能评价对比Table 2 Comparison of statistical parameters of different linear regression models (GAM,GWR and MGWR)
MGWR 模型的结果能直接表现出不同海洋环境因子的差异化作用尺度(表3),而GWR 模型仅能反映各变量作用尺度的平均状态。经典GWR 模型的带宽为65,为样本总数量的5.29%。MGWR 模型显示各海洋环境因子的作用尺度有很大差异,最佳拟合带宽由小到大依次为NPP、U55、S100、SST、V55。截距表示在其他环境因子确定的情况下,渔获率所受到的影响,其作用尺度即带宽为48,占样本总数量的3.91%,较低于其他变量的作用尺度。SST 的作用尺度也较小,为54,占样本总数的4.40%,说明海表面温度存在的空间异质性较大。S100 和U55 作用尺度均为48,占样本总数3.91%,表明这两个环境因子也具有较大的空间异质性。NPP 作用尺度最小,占样本总数量的3.51%,表明空间异质性极强。V55 带宽为1 227,属于全局尺度,其系数在空间上表现为平稳缓和,几乎不存在空间异质性。
表3 MGWR 与 GWR 带宽对比结果Table 3 Bandwidth comparison between classical MGWR and GWR models
在CPUE 与环境影响因素关系的实际讨论中,GWR 和MGWR 模型都考虑了空间非平稳性特点及空间尺度问题,但MGWR 模型为每个影响因素纳入了不同的空间尺度,充分利用每个变量的不同带宽(表3),对中西太平洋鲣的渔获率进行了更精确的回归分析。MGWR 模型的结果中,每个海洋环境因子有其特定的回归系数。各环境因子的具体回归系数信息见表4。
表4 基于MGWR的各环境因子的局部系数统计描述Table 4 Statistical description of MGWR local coefficient
本文中MGWR 模型各环境因子局部回归系数和显著性的空间分布(图1)显示,SST、S100、U55、V55 和NPP 5 种环境因子存在显著的空间差异。MGWR 各因子局部回归系数统计描述见表4。可以发现,SST 在大约2/3 区域对鲣渔获率有正向影响,1/3 部分为负向影响,170°E 为界东西两侧分别为负向影响和正向影响(图1-b),其中在赤道150°E 附近回归系数值最大。影响的平均值为0.162,且标准差较小,说明海表面温度每增加1 ℃,渔获率平均增加0.162 个单位。由显著性分布可知,170°E 以西海域,SST 对鲣渔获率正向影响的显著性较强,180°经线处负向影响最为显著。
图1 MGWR 模型各环境因子局部回归系数和显著性的空间分布Fig.1 Spatial distribution of local regression coefficients and significance of each environmental factor in MGWR model
S100 正向影响区域近73%,平均值为0.574,远大于其他因子,同时标准差也较大,因此变异系数相对较小,表明影响的空间差异较小。从系数空间分布图(图1-c) 可见,总体上S100 对鲣渔获率影响的显著性由南至北逐渐增强,其中有两个较为明显的区域,其一为赤道175°E 附近,另一显著区域位于赤道西太平洋“暖池”区域;负向影响主要分布在10°S 附近。
U55 局部回归系数取值介于-2.215~0.997,且变异系数较大,表明U55 对渔获率影响的异质性较强,负值比率约68%,因此主要为负向影响。总体上U55 对鲣渔获率的影响呈现环形分布状态,且负向影响范围更加广泛。标准差较大为0.504,说明局部回归系数差异也较大。由图1-d 可知,赤道170°E 附近为负向影响中心,影响程度由此向四周递减;在140°E—155°E、0°—10°S 内该系数变为正值,且根据显著性P值分布可知,此处U55 对鲣渔获率的正向影响最为明显,但影响范围较小。
V55 局部回归系数介于-0.141~0.553,平均值为0.069,标准差为0.109,正值比率超过80%,可见V55 在全局上对渔获率产生正向影响,由图1-e可知,V55 的负向影响区域较为分散;但从系数绝对值来看,其影响程度相对最小。
NPP 主要在65%研究区域内对渔获率产生正向影响,局部回归系数介于-0.787~1.806,平均值为0.172,标准差为0.351。如图1-f,170°E 两侧差异明显,在赤道150°E—160°E 附近NPP 显著性最显著,表明该处与鲣渔获率有较大相关性。170°E以东区域,局部回归系数主要为负值,对渔获率产生负向影响,且显著性较弱。
总体上,各环境因子对鲣渔获率影响的空间异质性程度可由变异系数表示,本研究中空间异质性程度依次为:U55>SST>NPP>S100>V55。由局部回归系数的正负值比例可知,V55 和S100 正负值相差比例及正值比例均较大,U55 的负值比例最大。各环境因子对鲣渔获率影响的空间异质性呈现东西向差异,SST 和NPP 对渔获率的影响主要受经度位置影响,总体上170°E 两侧会呈现不同的分布规律,S100 对渔获率的影响南北纬差异较为突出(图1)。同时,局部回归系数的绝对值大小表示各环境因子对鲣渔获率影响程度大小。本研究显示,各环境因子对渔获率的影响程度依次为:S100>U55>NPP>SST>V55。
由GAM、GWR 和MGWR 模型预测的中西太平洋鲣渔获率标准化后的空间分布图(图2)可见,不同模型的预测范围有所差异,其中GAM 模型预测的渔获率高值区介于140°E—178°W、5°S—2°N 之间,预测的渔获率低值区域主要分布在研究海域外围 (图2-a)。MGWR 模型预测的渔获率高值区见图2-b,分布在A1 和A2 海域范围内,渔获率低值分布范围主要在B1、B2、B3 及B4 附近(图2-b)。与GAM 模型相比,MGWR 模型预测的渔获率高值区范围较小,而低值区范围相对较广。GWR 模型预测的渔获率高值分布区域与MGWR 模型的结果几乎一致,渔获率低值区主要为C1、C2 和C3(图2-c)。与GWR 模型相比,MGWR 模型预测的渔获率高值区域范围更集中。总体而言,由3 个模型预测的渔获率空间分布与2005—2019 年CPUE空间分布(图2-d)可知,与GAM 模型相比,MGWR 和GWR 模型能更好地反映渔获率的分布状况。此外,由于各海洋环境因子具有空间异质性特点,GAM 模型无法明确体现其对结果产生的影响,而MGWR 和GWR模型均能较好地体现出这一特点。
鲣是大洋性中上层高度洄游性鱼类[24],其资源丰度、集群、资源的分布及洄游与海洋环境息息相关[25]。本研究中MGWR 模型结果显示(图1),SST、S100、U55、V55、NPP 5 种环境因子对鲣渔获率的影响均存在显著的空间差异。S100 对鲣渔获率主要为正向影响,S100 影响下出现两个较为明显的正向区域,主要原因为赤道170°E 附近受热带东太平洋信风影响,生成了巨大的涌升流,因此形成了低温、高盐的冷舌区域[26];S100 影响下另一个高值区域和赤道西太平洋“暖池”区域重合度高,“暖池”区域有高温、低盐的特征,当盐度增加时,净初级生产力提升,从而为鲣提供食物资源[27],导致产生更显著的空间聚集效果(图1-c)。
U55 对研究区域内鲣渔获率主要为负向影响(图1-d)。根据洋流分布可知,本研究区域有多条东西向暖流经过,大部分区域都受到了影响,处于洋流流经区域,因此在区域内东西向海流速度对渔获率分布差异产生了较大影响。本研究中U55 正向影响明显的区域处于赤道附近,受赤道暖流影响,洋流方向引导鱼群聚集于此[28]。
NPP 对渔获率主要为正向影响。浮游动植物的生长通常会受到NPP 的影响[29],NPP 值越大的区域,其环境条件对鲣的生存越有利,不同季度冷暖水团的锋面位置会发生变动,NPP 通常会随冷暖水团的锋面移动[30],因此本研究中NPP 主要对鲣渔获率存在正向影响。
SST 在170°E 以西对渔获率有显著的正向影响。已有研究显示,中西太平洋热带海域是一个“暖池-冷舌海洋生态系统”海域[21,31],在图1-b中,正向影响最为显著的区域基本在中西太平洋的“暖池”区域内,空间分布范围较为广泛,整体资源热度均较好,其整体环境特征是高温和相对较低的叶绿素浓度,并且“暖池”海域靠近陆地及岛屿,因此近岸上升流也提供了较充足的初级生产力[32]。同时在“冷舌”附近渔获率也受到了一定程度的正向影响,主要是受中东太平洋赤道上升流影响,锋面区域也产生了较丰富的初级生产力。金枪鱼鱼群往往在锋面地带获得丰富的食物资源[29,31],这和本研究的结果基本吻合。
总体上,V55 对鲣渔获率影响的空间异质性程度相对最小(表4),空间异质性的影响范围较小,在正负影响的差异变化上较为平稳缓和(图1-e);但V55 对渔获率的正向影响相对最为明显,主要体现在靠近巴布亚新几内亚的海域及180°附近(图1-e)。U55 空间异质性程度最大,但目前较少有研究使用海水流速等动力学因子进行分析,未来的研究建议考虑相关的海洋动力学因子。SST 和NPP 也有较明显的影响。S100 对渔获率的影响在纬度变化上也有东西向差异,但南北纬向差异更为突出,NPP和SST 对渔获率的影响主要受经度位置影响,总体上170°E 两侧会呈现不同的分布规律。这是由于170°E 以东海域受热带东太平洋信风影响,从而产生了涌升流,形成了温度低、盐度高、初级生产力高的冷舌区。暖池边缘生成的由底层海水上泛且带有大量营养物质的辐合区,聚集了大量的微型浮游动物及浮游植物,使鲣资源拥有了丰富的饵料[33],形成了良好的产卵场和索饵场。在170°E 以西鲣的空间聚集区附近海域,海表面温度相对较高,西部暖锋有两个“暖池”向东凸出,有利于锋面的形成,产生的丰富食物资源为鲣提供了良好的生存条件,鲣的资源量因此较为丰富。
本研究中MGWR 模型各环境因子回归系数的绝对值显示,各环境因子对渔获率影响程度大小不同。5 个影响因子中,S100 对渔获率的影响程度最大,其次为U55、NPP 和SST,V55 的影响程度相对较小。盐度对鲣渔场影响的研究表明,盐度在海洋变化中充当重要的物理参数[34-35]。大洋环流中海水盐度是有效的示踪物,且通过对海水密度的调节变化,在海洋热量的垂直输送上产生了重要影响[36],海表面温度异常也间接受到影响[37]。本研究中热带中西太平洋海域水分和温度条件充足,盐分对表层营养物质至关重要,从而对鲣资源产生重要影响。同时混合层深度受温度和盐度作用,对鲣资源也产生了一定影响[38]。U55 对渔获率大小的影响程度也相对较高。Meehl[28]认为,东西向海水流速在热带地区有最大的季节性变化,海流的变化会引起盐度、溶解氧、温度等因素发生改变,使诱饵产生聚集,海流方向可以引导鱼类游向诱饵来源,鱼群更容易定位诱饵。受信风的影响,东太平洋产生了强烈的涌升流,形成了盐度高、温度低等特征的冷舌区域。冷暖流相遇形成的交汇海域浮游动植物较为充足,是鲣理想的索饵场。因此,本研究认为NPP 也是影响鲣渔获率的重要因子之一。张小龙等[39]认为,由于148°E—165°E、6°N—6°S 附近海域存在两个向东凸出的“暖池”易形成温度锋面,为鲣的生存提供了丰富的饵料资源,形成了明显的空间聚集现象,因此SST 在一定程度上也影响了中西太平洋鲣渔获率。以往研究认为SST 是影响鲣渔获率的重要因子,但本研究采用平均状态下的环境和资源数据,因此SST 的影响程度较次于S100 和NPP。同时从空间形态观察来看,S100、NPP 和SST 这三者区域差异的影响模式相似。
本研究中,由于MGWR 模型考虑了环境因子的多尺度效应,结果发现NPP 作用尺度较小,其次为S100 和U55,SST 的作用尺度较大,V55 最大,为全局尺度。NPP 是主要表现食物的因子,其作用尺度较小可能是由于鲣资源捕食范围受该区域冷暖水团的影响,而冷暖水团边界范围较狭窄,所以其影响的空间范围也相对较小。鲣是喜温性鱼类,主要分布在水温高于29.5 ℃的区域内[34],本研究区域内暖水团范围较广,因此SST 的作用尺度较大。近岸上升流和盐度锋面均对营养盐有一定影响[32],但由于位置较为固定,因此S100 的作用尺度也受到了一定的空间范围限制。
总体上,SST、NPP 对鲣渔获率空间分布的影响相似,受到经度位置的影响,170°E 两侧呈现不同的规律,其中170°E 以西主要为正向影响,且具有明显的岛屿属性;170°E 以东主要为负向影响。这是由于170°E 以东海域在热带东太平洋信风的作用下,形成温度低、盐度高的冷舌区域;170°E 以西海域,正向影响最为显著区域基本在中西太平洋的“暖池”区域内,整体资源热度均较好,整体上温度较高、叶绿素浓度相对较低,并且该区域靠近陆地和岛屿,因此近岸上升流也为其提供了丰富的营养物质。S100 南北向差异较为明显,而V55 没有明显的东西或南北向差异。U55 正向影响区域较小,主要集中在热带岛屿附近,原因可能是海流的变化引起岛屿附近盐度、溶解氧、温度等因素发生变化,使饵料产生了聚集[28]。
目前,有关环境与鲣关系的研究主要采用MaxEnt 模型[12]、GAM[40]模型等,但这些模型均未考虑空间异质性问题。GWR 模型可以通过结合不同地理位置的空间属性差别,来探索不同环境因子在时空上的差异性影响,并获得海洋环境与资源间的时空异质性特征。本研究在GWR 模型基础上引入MGWR 模型,加入了多尺度概念,能捕捉到不同变量的不同影响尺度,考虑了不同环境因子在影响尺度上的差异性,从而避免了捕获太多的噪声和偏误。因此,是否考虑影响因素的空间尺度会对模型的结果和分析产生巨大影响。由于本研究采用平均状态的环境和资源数据,所用的MGWR 模型未顾及时间序列对海洋资源的影响,建议未来的研究在此基础上加入时空地理加权回归模型,综合时间和空间两个维度进行分析。