环境因子对东南太平洋智利竹䇲鱼渔场时空分布异质性影响

白思琦，邹晓荣, ，张 鹏，丁 鹏

（1.上海海洋大学海洋科学学院，上海 201306；2.国家远洋渔业工程技术研究中心/大洋渔业资源可持续开发省部共建教育部重点实验室/农业农村部大洋渔业资源环境科学观测实验站，上海 201306）

智利竹䇲鱼 (Trachurus murphyi) 是一种具有较高商业价值的高度洄游性经济鱼种[1-4]，主要分布在东南太平洋的中上层海域，包括厄瓜多尔和秘鲁等沿海国水域以及智利专属经济区向西延伸的部分外海水域[5-8]。其产量在相关沿岸国家及一些远洋渔业国家的捕捞产业中具有举足轻重的地位[9-12]。智利竹䇲鱼渔业自20 纪70 年代被大规模开发以来，产量多年占据世界前列[13]，1995 年达到高峰 (4.96×107t)。近20 年来，智利竹䇲鱼资源量由于全球气候长期变化和人类捕捞强度的增大而不断降低，渔场的年际变动愈加显著[14]。

渔场资源变动与海洋环境息息相关[15-21]，对此国内外学者利用不同模型研究海洋环境因子对智利竹䇲鱼渔场时空分布的影响，Laura 等[16]通过神经网络模型、汪金涛等[22]利用主成分和BP 神经网络方法、晋伟红[23]利用偏最小二乘法计算模型、陈春光[24]和Li 等[25]通过建立不同权重的栖息地适应性指数模型、杨香帅等[14]采用灰色关联度分析方法研究了智利竹䇲鱼渔场时空变动与海洋环境因子间的关系，这些模型的全局性研究会忽略资源分布与环境因子的特殊区域。

海洋生态系统具有高度复杂性，鱼类资源由于受到多种海洋环境因子的交互作用产生动态响应机制，出现空间上的变动，且该变动具有空间异质性[26]。目前对智利竹䇲鱼资源的研究均偏向于静态，而地理权重回归模型 (Geographically weighted regression,GWR) 将采样点间的空间距离考虑在内，对每个空间位置拟合一个局部回归方程，得出各环境因子与单位捕捞努力量渔获量 (CPUE) 关系的局部回归系数，便于发现资源分布的“热点区域”[27-29]。Windle 等[28]首次利用GWR 模型分析渔业数据，发现环境因子对渔业资源的空间分布有异质性影响，目前该模型已广泛应用于陆上资源空间分布的研究，在渔业资源空间分布研究中也有较为广阔的应用前景。本研究运用GWR 模型对影响智利竹䇲鱼种群生物资源丰度的环境因子进行对比分析，以期为我国智利竹䇲鱼渔船的生产提供参考。

1 材料与方法

1.1 数据来源

智利竹䇲鱼捕捞活动数据由中国远洋渔业协会大型加工拖网项目技术组提供，数据来源于大型拖网渔船“开富号”和“开裕号”，两船为相同类型的中层拖网船，使用的渔具规格为1 632 ×444.00 m(277.40 m)，网板面积为12 m2。船长、型深、船宽、总吨位、主机功率分别为91.10 m、9.20 m、20 m、7 671 t、5 920 kW。捕捞日志记录数据包括作业日期、作业水深、起放网时间、作业位置以及渔获量等。作业区域为73°W—92°W、19°S—46°S。由于每年9 月到翌年3 月为生产淡季，出于对数据充足性的考虑，时间范围选择为2012—2018 年的4—8 月。去除渔捞日志中捕捞数据记录有误的数据，有效数据为1 880 个。环境数据选取海面温度 (Sea surface temperature,SST)、叶绿素a浓度 (Chlorophyllaconcentration,Chl-a)、海面温度异常值 (Sea surface temperature anomaly,SSTA)和海面高度异常值 (Sea surface height anomaly,SSHA)及作业水深，海面温度异常值是指实际海面温度按照时间序列减去该时间序列的平均海面温度所得到的海面温度差值，海面高度异常值是指实际海面高度按照时间序列减去该时间序列的平均海面高度所得到的海面高度差值。SST、Chl-a、SSTA 和SSHA 来源于美国国家大气局 (NOAA)Oceanwatch 网站 (www.oceanwatch.pifsc.noaa.gov/)，作业水深为实测数据，时间范围为2012—2018 年4—8 月，空间分辨率为0.125°×0.125°，时间分辨率为周。

1.2 数据处理及分析

1.2.1 CPUE 计算 CPUE (t·h−1) 是用于表示智利竹䇲鱼资源密度的主要指标[14]，计算公式为：

1.2.2 渔场重心计算 分别统计各个年度每一经度和纬度的渔获量，对各周作业的渔场重心进行计算，并分析比较渔场空间分布的变化，以0.5°×0.5°的空间分辨率按周进行统计，作业渔场的经纬度重心计算公式[30]为：

式中X表示渔场经度重心，Xi表示第i渔区的经度中心，Ci表示第i渔区的渔获量，Ei表示第i渔区的捕捞努力量；Y表示渔场纬度重心，Yj表示第j渔区的纬度中心，Cj表示第j渔区的渔获量，Ej表示第j渔区的捕捞努力量。

1.2.3 空间自相关分析 用于表示空间自相关性方法众多，Moran's I 指数为国际中最广泛使用的方法，该指数用于衡量空间要素的相互关系，数值范围为 [−1,1]，大于0 表示呈正相关，反之则为负相关，经检验等于0 则表明不存在空间相关性[31-33]。该方法表达式为：

1.2.4 地理权重回归模型 GWR 模型是将地理学第一定律深度融入局部空间统计方法的模型，通过对独立抽样的捕捞站点分别进行模型的回归分析计算，得到与空间位置相对应的空间回归系数，随着空间位置不同而变化的参数估计量化表征空间关系异质性特征[29]。若智利竹䇲鱼CPUE 与环境因子之间的关系存在空间异质性，则GWR 模型中变量的关系可能存在空间位置函数，因此采用反映局部空间分异的GWR 模型能在不同局部空间进行参数估计[34]。GWR 模型是全局线性回归模型的扩展，根据每个空间位置周围一定范围内的作业点，都可得到一个局部回归模型的结果数据[27-28,35-36]。与传统线性回归模型的区别在于GWR 模型中用于拟合局部回归模型各采样点的权重值不是假定相等的，而是与距离回归点的距离成反比，采样点权重值的确定常常采用高斯和双平方权重函数[35]。GWR模型公式[34]表示为：

式中yi表示智利竹䇲鱼CPUE，βj(j=0、1、2、3、4、5) 表示回归系数；(ui,vi) 表示地理空间中第i个回归点的空间坐标；ESSTi表示地理空间中第i个回归点的海面温度；EChlai表示地理空间中第i个回归点的Chl-a浓度；ESSTAi表示地理空间中第i个回归点的海面温度异常值；ESSHAi表示地理空间中第i个回归点的海面高度异常值；Edepthi表示地理空间中第i个回归点的作业水深。

1.2.5 显著性检验 对地理权重回归模型进行显著性检验和模型误差项空间自相关的显著性检验，常用的检验方法为赤池信息准则 (Akaike information criterion,AIC) 检验[37]、蒙特卡罗检验[38]和F检验[39-40]。本文选取AIC 检验，AIC检验多用于比较不同模型的显著性，当两个模型AIC 值相差超过3 时，较低AIC 值的模型具有更好的拟合效果。表达式为：

式中n是样本数，σ是误差项估计的标准离差，tr(S)是地理加权回归模型的投影矩阵S 轨迹，是带宽的函数。

本研究使用ArcGIS 10.6 软件来进行数据空间自相关性分析、GWR 模型建立及模型结果的显著性检验。

2 结果

2.1 渔场重心时空变化

2.1.1 CPUE 及渔场重心月际变化 智利竹䇲鱼渔场2012—2018 年CPUE 月际变动见图1，智利竹䇲鱼渔场2012—2018 年4—8 月CPUE 总体基本呈现为逐月下降的趋势，仅2012、2013、2014年基本呈上升趋势；8 月CPUE 波动最小，2015 年数值最大 (10.64 t·h−1)，2017 年数值最小 (6.11 t·h−1)；4 月CPUE 波动最大，2017 年数值最大 (14.08 t·h−1)，2013 年数值最小 (3.76 t·h−1)。智利竹䇲鱼渔场2012—2018 年月际变动见图2，渔场逐月向北移动，4 月渔场主要分布在78°W—87°W、40°S—45°S；5 月主要分布在78°W—90°W、40°S—45°S；6 月主要分布在78°W—88°W、38°S—44°S；7 月主要分布在76°W—90°W、33°S—41°S；8 月主要分布在73°W—88°W、19°S—39°S，且开始分为两个区域，一部分逐渐向西北偏移，一部分向东北偏移；8 月渔场位置达到最北端，5 月达到最南端。

图1 2012—2018年月平均单位捕捞努力量渔获量变化图Figure 1 Monthly CPUE change from 2012 to 2018

图2 渔场重心月际变化Figure 2 Monthly gravity change in fishing ground

图3 渔场重心年际变化Figure 3 Anually gravity change in fishing grounds

2.1.2 CPUE 及渔场重心年际变化 智利竹䇲鱼渔场2012—2018 年CPUE 年际变动见图1，2012—2018 年年际CPUE 基本呈逐年上升的趋势，2012 年CPUE 值最小(4.82 t·h−1)，2015 年CPUE 值最高(10.94 t·h−1)。智利竹䇲鱼渔场2012—2018 年年际变动见图3，2012 年渔场主要分布在75°W—88°W、27°S—45°S；2013 年主要分布在77°W—81°W、36°S—45°S；2014 年主要分布在75°W—82°W、27°S—44°S；2015 年主要分布在77°W—90°W、34°S—45°S；2016 年主要分布在74—86°W、28°S—44°S；2017 年主要分布在73°W—84°W、20°S—45°S；2018 年主要分布在73°W—90°W、22°S—44°S；2017 年渔场位置达到最北端，2015 年达到最南端，2015 年渔场位置最为集中。

2.2 空间自相关性分析结果

Moran's I 指数结果显示CPUE 与SST、Chl-a、SSHA、SSTA、作业水深 (Depth) 均具有显著的正空间自相关性 (Moran's I 指数>0，表1)，且该变量空间分布格局较为清晰，可以进行GWR 模型的建立。

表1 不同环境因子地理权重回归模型回归效果的比较Table 1 Comparison of regression effects of GWR models with different environmental factors

2.3 CPUE 与环境因子关系的空间异质性

GWR 模型结果显示智利竹䇲鱼CPUE 与SST、Chl-a、SSTA、SSHA 和作业水深在空间分布上具有异质性 (图4)。回归系数值的大小表明GWR 模型中该环境因子对渔场时空分布影响的强弱，正值表明该环境因子对CPUE 存在正向影响，负值表明存在反向控制作用。SST、Chl-a、SSTA、SSHA、作业水深的局部回归系数值都存在正负符号的改变，即不同区域环境因子对CPUE 的影响不同，数值的大小表明了这5 个环境因子在空间上均对CPUE 的影响显著；整体来看，基于GWR 模型的CPUE 与SSTA 回归系数比Chl-a、SSTA、SSHA、作业水深大，表明SSTA 对于东南太平洋智利竹䇲鱼渔场分布的影响大于其他4 个因子；CPUE与SST、Chl-a、SSTA、SSHA、作业水深的回归系数大部分为正值，即SST、Chl-a、SSTA、SSHA、作业水深对CPUE 均具有正向影响。

图4 地理权重回归模型评估的单位捕捞努力量渔获量与海面温度、叶绿素 a 浓度、海面温度异常值、海面高度异常值和作业水深关系的局部回归系数值Figure 4 Local regression coefficients of SST,Chl-a,SSTA,SSHA and gear depth by GWR model

自然环境条件具有显著的空间差异，因此CPUE 与SST、Chl-a、SSTA、SSHA、作业水深也是呈现显著的空间异质性。拟合优度 (R2) 和局部拟合优度 (校正的R2) 越大则表示对变量的解释率越大，校正的R2计算是按照分子和分母的自由度对它们进行正规化，具有对模型中变量数进行补偿的效果。SSTA 基于GWR 模型回归的R2值最大，为0.74，SSHA 的校正R2值最大，为0.59，SSTA 的校正的R2值次之，为0.58 (表1)。结果表明SSTA 对东南太平洋智利竹䇲鱼渔场时空分布影响最大。

2.4 AIC 检验结果比较

AIC 检验见表1，结果显示，5 种环境因子中SST 基于GWR 模型回归的AIC 值最小 (1 022.08)，虽然其拟合优度较低，但拟合精度最好。

3 讨论

3.1 东南太平洋智利竹䇲鱼渔场时空变动

本研究发现2012—2018 年4—8 月智利竹䇲鱼渔场分布在73°W—90°W、20°S—45°S，2017年渔场位置达到最北端，2015 年达到最南端且渔场位置最为集中，为77°W—90°W、34°S—45°S。魏广恩等[41]研究认为渔场位置可以用渔场重心来判断。西风漂流带为45°S—60°S，邹晓荣[42]研究认为，西风漂流带为智利竹䇲鱼南部的渔场屏障，本文研究结果与其一致。在渔场月际变化趋势上，渔场逐月向北移动，与秘鲁寒流从南向北的流向有关，8 月渔场分为两部分，一部分逐渐向西北偏移，一部分向东北偏移，该变动可能与秘鲁寒流在该区域分成两部分，一部分流入南赤道暖流，一部分流向东澳大利亚暖流有关。2012—2018 年4—8 月CPUE 基本呈逐月下降的趋势，8 月CPUE 最小，4 月CPUE 最高，表明在生产旺季中4 月为生产最高，可为今后东南太平洋智利竹䇲鱼捕捞生产提供科学的管理依据。

21 世纪初，由于环境变化和捕捞导致资源骤减，2006 年成立了南太平洋区域渔业管理组织，限制各国的捕捞水平，旨在恢复智利竹䇲鱼资源量，在时间年际变化趋势上，智利竹䇲鱼渔场2012—2018 年CPUE 基本呈逐年上升的趋势，说明了南太平洋区域渔业管理组织总捕捞许可量制度管理初见成效。有研究认为厄尔尼诺和拉尼娜现象会引发环境发生异常改变，进而导致海洋环境生产力的变化，资源分布产生差异[14]，研究表明2015 年为厄尔尼诺年[12]，受其影响，CPUE 值最高，且探捕报告调查表明生产时间和产量均增加，可能与厄尔尼诺现象有关。

3.2 东南太平洋智利竹䇲鱼CPUE 与环境因子关系

本文研究了2012—2018 年4—8 月环境因子对东南太平洋智利竹䇲鱼渔场时空分布的异质性影响，结果表明，环境因子中SST、Chl-a、SSHA、SSTA、作业水深存在显著的时空差异性，表明了智利竹䇲鱼渔场时空分布与该5 种环境因子存在异质性关系，由于复杂的海洋环境以及较大的研究范围与海洋环境因子的相互作用导致CPUE 与环境因子存在空间异质性关系。结果表明SSTA 对于东南太平洋智利竹䇲鱼渔场时空分布影响最大，说明厄尔尼诺和拉尼娜现象更能影响渔场的时空分布；CPUE 与SST、Chl-a、SSTA、SSHA、作业水深的回归系数多为正值，即该5 种环境因子对CPUE具有正向影响，在78°W—85°W、35°S—45°S 正向关系尤其显著，从整体来看，基于GWR 模型回归系数的对比，海面温度和SSTA 对渔场时空变动的影响大于其他环境因子，相关空间格局的结果表明，渔场从南到北，自东向西移动。牛明香等[9]基于GAM 模型将各个环境因子综合起来定量评估其重要程度及影响，但GAM 模型得出的全局统一规律可能会忽略资源分布与环境因子的特殊区域，GWR 模型的使用，考虑了数据空间关系的局部特征，可研究渔场的动态分布，便于发现智利竹䇲鱼资源空间分布的“热点”海域，目前GWR 模型多用于陆地资源的研究，未来其在海洋研究上具有广阔的应用前景，也能对未来渔业资源的开发和利用提供参考依据。

本研究仅分析了部分环境因子，还有许多环境因子例如海水盐度、溶解氧、海流和食物来源等，均对CPUE 产生一定影响，之后可将其融入到模型中，增加模型精确度，得到更为客观的分析结果；本文数据为我国大型拖网渔船的商业捕捞数据，缺少他国渔船捕捞数据，在今后的研究中需要对结果进行进一步验证；GWR 模型不适合预测调查海域以外的资源分布；在经济、自然资源管理、医学等多个领域的研究趋向于平面型，在地理上仅部分研究可呈立体型，例如不同海拔具有不同温度、降水、土壤环境等，与海洋环境相似，但部分海洋生物具有垂直移动特性，使用该模型时需要考虑具体情况或改变相关影响因子。