毛里塔尼亚海域短线竹䇲鱼时空分布与海洋环境的关系

宋利明，徐双泉，许 回，陈明锐，Ebango Ngando Narcisse

（1.上海海洋大学海洋科学学院，上海 201306；2.国家远洋渔业工程技术研究中心，上海 201306；3.大洋渔业资源可持续开发省部共建教育部重点实验室，上海 201306；4.远洋渔业协同创新中心，上海 201306）

短线竹䇲鱼（Trachurus trecae）为中底层鱼类，一般活动在20～100m的近海底水域，有时也在中上层甚至近表层出现，适宜水温为15～22℃，广泛分布在东大西洋的摩洛哥至安哥拉海域［1］。毛里塔尼亚西邻大西洋，处在几内亚湾暖流和加纳利寒流的交汇地带，海域初级生产力比较高，拥有丰富的渔业资源［2］。FAO的统计结果显示［3-4］，2007年毛里塔尼亚海域主要的小型中上层鱼类中，短线竹䇲鱼占总渔获量的比重最多，为31%；2013年中东大西洋的短线竹䇲鱼产量主要集中在毛里塔尼亚沿海，约占82%。因此对毛里塔尼亚海域短线竹䇲鱼的研究十分重要。目前国内学者如牛明香等［5-6］研究了智利竹䇲鱼（Trachurusmurphyi）的渔场与时空变动的关系，李显森等［7］、化成君等［8］和汪金涛等［9］分别对智利竹䇲鱼的资源分布特征、时空分布与环境因子的关系、渔情预报技术作了研究。RUMOLO等［10］通过同位素和胃含物分析，研究了西西里海峡大西洋竹䇲鱼（Trachurus trachurus）的空间分布与摄食行为的关系，GUSHCHIN和CORTEN［11］研究了毛里塔尼亚海域大西洋竹䇲鱼和短线竹䇲鱼的摄食特性，ARKHIPOV［12］研究了毛里塔尼亚海域鱼类浮游生物的季节和年际变化，班锐［13］利用广义可加模型探讨了时空和环境因子对毛里塔尼亚海域竹䇲鱼渔场的影响，但毛里塔尼亚海域的竹䇲鱼包括了短线竹䇲鱼在内的多种竹䇲鱼类，其具体的生活方式和时空分布特征存在差异［1］，且其研究中使用的环境数据来源于卫星遥感，分辨率较低。分位数回归模型最早由KOENKER和BASSETT［14］提出，已在研究鱼类的栖息地综合指数方面得到了应用［15-16］。本文根据2017年在毛里塔尼亚海域测定的海洋环境数据以及短线竹䇲鱼的渔获率统计资料，用分位数回归的方法探究短线竹䇲鱼时空分布与海洋环境的关系，旨在进一步了解短线竹䇲鱼在毛里塔尼亚海域的时空分布特征，为提高我国渔船在此海域的作业效率及毛里塔尼亚对短线竹䇲鱼等小型中上层渔业资源养护与管理提供参考。

1 材料与方法

1.1 调查方法

1.1.1 调查船及基本信息

“福远渔097”灯光围网船总长49.95 m，型宽8.60 m，型深4.00 m。设计吃水3.30 m，设计排水量896.90 t，主机功率为828 kW，航速为6～11 kn，鱼舱容量190m3，油舱容积120m3，淡水舱容积47m3。航区为国际I区。

1.1.2 调查时间、海域及站点分布

调查时间为2017年9月20日至12月31日。调查海域为：16°28′～17°59′W、17°50′～20°58′N。共119个调查站点（部分站点位置非常接近，图中重叠在一起），随机选择75%的（90个）站点的数据用作建模，25%的（29个）站点数据用以验证模型，站点分布见图1。部分建模站点和验证站点出现重合，原因是这些站点的下网位置几乎相同。

1.1.3 调查的渔具渔法

“福远渔097”所用围网网具主体由网衣、包铅绳、浮子、上纲和下纲组成，网口网衣拉直周长约1 643.4 m，网口至囊端拉直长度约116 m，上纲长度约809.84 m，下纲长度约761.44 m。

1.1.4 调查仪器、方法及内容

利用渔船上的全球定位仪记录每天放网开始时的经纬度，同时记录渔获重量。调查站点预先确定，但实际调查位置与计划站点存在一定偏差。每日放网后或起网前，使用加拿大型号为XR-620的多功能水质仪，测量调查站点海表面的叶绿素a浓度、海表面温度以及盐度数据。

1.2 数据处理方法

1.2.1 CPUE的计算方法

单位捕捞努力量渔获量（catch per unit effort，简称CPUE）是用来表示鱼类资源丰度的相对指数，是对某个渔业资源进行评估的基本内容之一［17］。船上所得到的直接数据为渔获量（t），建立关系模型前，需要计算短线竹䇲鱼的CPUE值，计算方法如下：

式（1）中，ci为在i站点的渔获量（t），ni为i站点的网次数量。

1.2.2 海洋环境因子值的计算方法

海洋环境因子数据为各个站点10 m以浅所测得的各海洋环境因子数据的算术平均值［15］。

1.3 建立短线竹䇲鱼CPUE预测模型

分位数回归建模的具体方法可以参照文献［15］。分位数回归分析采用由美国地理调查局中陆生态科学研究中心开发的统计学软件Blossom。

根据90个建模站点的CPUE、叶绿素a浓度（Ch）、海表面温度（T）、海表面盐度（S）及其交互作用项建立短线竹䇲鱼CPUE与海洋环境因子的关系模型，其一般形式为：

图1 调查站点分布图Fig.1 Distribution map ofmodeling sites and validation sites

式（2）中，C为常数项，a-f为自变量的系数，ChT表示叶绿素a浓度和海表面温度的交互作用项，ChS表示叶绿素a浓度和海表面盐度的交互作用项，TS表示海表面温度和海表面盐度的交互作用项，ε为误差项。i表示站点号数，如Chi表示在i站点的叶绿素a浓度。

1.4 栖息地综合指数（IHI）计算方法

根据关系模型，得出建模站点和验证站点的预测CPUE，按照下列公式计算得出短线竹䇲鱼在不同站点的IHI值：

式（3）中，CPUEij是指j组（j=1为建模站点；j=2为验证站点）i站点的CPUE预测值，CPUEmaxj是指j组中CPUE预测值中的最大值。

1.5 IHI模型预测能力的验证和评价

将建模站点和验证站点的环境数据值分别代入模型中，使用SPSS 23.0软件，通过Wilcoxon（符号秩）对建模站点和验证站点的CPUE实际值与预测值进行检验［18］，根据P值的大小判断两者间是否存在显著性差异，以此评估模型的预测能力。通过计算短线竹䇲鱼验证站点的CPUE实际值与预测值之间的Spearman相关系数（秩相关系数）大小，结合双尾检验对IHI模型的预测能力进行评价。

用Arcmap10.6作图得出建模站点和验证站点的IHI等值线与实测CPUE的分布图。将29个验证站点的实测CPUE值除以其中的最大值得出IHI′［式（4）］，用Excel作出IHI′和IHI的变化趋势图，定性评价模型的预测能力。

式（4）中，CPUEi2是指验证站点i的实测CPUE，CPUEmax2是指验证站点实测CPUE中的最大值。

2 结果与分析

2.1 分位数回归得出的预测模型

利用Blossom建立CPUE与海洋环境的关系模型，将每次拟合得出的模型系数和对应的P值记录下来，取分位数（用Q表示）最高的模型为最佳模型。结果如表1所示。

拟合得出的最佳预测模型如下：

由式（5）可知，影响短线竹䇲鱼CPUE的海洋环境因素包括海表面温度、叶绿素a浓度、海表面温度与叶绿素a浓度的交互项，P值依次为0.028、0.004和0.032。

由自变量系数的大小可得，对CPUE影响最大的是海表面温度，其次是叶绿素a浓度，再次是海表面温度和叶绿素a浓度的交互作用项。

2.2 IHI模型验证

对建模站点和验证站点的CPUE实测值与预测值进行Wilcoxon（符号秩）检验，结果见表2，两者的P值均高于显著水平0.05，说明建模站点和验证站点的CPUE实测值与预测值之间均无显著性差异。

对29个验证站点的实测CPUE和预测CPUE的相关系数和显著性进行双尾检验，得出实测CPUE与预测CPUE的Spearman相关系数达到0.832，相伴概率小于0.01，这说明显著性水平为0.01时，实测CPUE和预测CPUE之间高度相关，模型可靠。

表1 不同分位数下各变量回归系数秩得分检验结果Tab.1 Rank score test for regression parameters under different quantiles

表2 建模站点和验证站点的实测CPUE与预测CPUE间W ilcoxon检验结果Tab.2 Result of W ilcoxon test between the observed and predicted CPUE ofmodeling sites and validation sites

2.3 IHI模型的定性评价

用Arcmap10.6作图得出的建模站点的IHI等值线与实测CPUE分布见图2。验证站点的IHI等值线与实测CPUE分布见图3。

本次调查期间，短线竹䇲鱼的CPUE在0～3 t·次-1，因此在作图和分析时，认为CPUE小于1 t·次-1时为CPUE的较低值，CPUE大于2 t·次-1时为CPUE的较高值。由图2可观察到，实测CPUE与IHI具有一定的相关性。IHI值在小于0.40时，实测CPUE值小于1 t·次-1；IHI值在0.40～0.56之间时，CPUE值随IHI值的增大而相应增大，大部分在1～2 t·次-1之间，大于2 t·次-1的情况也有出现，但分布零散，具有偶然性；IHI值大于0.60时，CPUE较高值出现的频率明显增加，大于2 t·次-1的情况集中出现。

观察图2可得，CPUE较低值集中分布在16°15′～16°50′W、16°50′～18°50′N，对应的IHI值低于0.56；CPUE较高值集中分布在17°20′～17°45′W、20°15′～20°50′N，对应的IHI值大于0.60。

图2 建模站点IHI等值线和实测CPUE分布Fig.2 IH I contour and CPUE distribution for modeling sites

图3中，验证站点的CPUE与IHI呈现出的相关性较为显著。当IHI值低于0.40时，对应的海域为16°15′～17°40′W、18°05′～20°05′N，CPUE较低值在这一区域零散分布，大于1 t·次-1的只有1个，其余都在1 t·次-1以下。当IHI值大于0.50时，对应的海域为17°20′～17°45′W、20°15′～20°50′N，CPUE较高值集中出现，CPUE全部大于1 t·次-1，大于2 t·次-1的验证站点也都在上述区域出现。这说明IHI对CPUE的预测具有良好的效果。

29个验证站点的IHI′和IHI的变化趋势见图4。可以观察到，IHI的变化趋势与IHI′高度一致，基本呈现正相关，侧面反映了预测CPUE的变化趋势与实测CPUE十分吻合。这可说明，建立的模型能较好地预测CPUE。

图3 验证站点IHI等值线和实测CPUE分布Fig.3 IHI contour and CPUE distribution for validation sites

图4 验证站点IH I′和IH I的变化趋势Fig.4 Trends of IHI′and IHI for validation sites

综上，IHI值较高的海域对应的实测CPUE较高，因此基于分位数回归拟合得出的关系模型能较好地反应CPUE与海洋环境因子的关系，预测能力良好。

3 讨论

3.1 用分位数回归模型及相关统计检验方法的依据

分位数回归可提供不同分位数下的回归方程，反映了自变量对因变量的某个特定分位数的边际分布，因此能清楚地阐释因变量的整个分配，甚至数据异质性的问题也能解决［19］。本文所使用的样本数据时间跨度只有3个多月，样本数据总量也只有119个，加上部分站点数据用作验证模型，实际上用于建模的数据并不多，因此样本数据并非呈正态分布，而是呈正偏态分布。当研究的数据集合呈现非正态分布时，选择分位数回归能够具体地反映因变量的平均水平［20］。利用SPSS对样本数据进行非参数检验，当样本数据呈现正态分布时，可用t检验，非正态分布或无法确定分布状态时，用Wilcoxon秩和检验［21］；进行相关性分析时，样本数据呈正态分布用Pearson检验，非正态分布使用Spearman检验［22］；未知呈正相关和负相关，使用双尾检验。

3.2 海洋环境因子的影响与时空分布

由基于分位数回归得出的最佳关系模型可以看出，海洋环境因子对短线竹䇲鱼CPUE的影响程度为：海表面温度＞叶绿素a浓度＞海表面温度与叶绿素a浓度的交互作用，海表面盐度对CPUE的影响没有在关系式中被表征出来。

渔场的形成与分布是由各种海洋环境因子的综合作用共同决定的。在所有非生物环境因子对鱼类资源分布的影响中，牛明香等［6］指出，水温是对渔场形成最为重要的因素，海表面温度的变化能够影响竹䇲鱼的分布和洄游等。本文的结论同样为海表面温度对短线竹䇲鱼的CPUE影响最大，这也符合实际情况：毛里塔尼亚海域处在几内亚湾暖流和加纳利寒流的交汇地带，有利于渔场形成，寒暖水团对应着温度的变化，因此海表面温度对短线竹䇲鱼的影响尤其显著。浮游生物是短线竹䇲鱼等中上层鱼类的食物来源，也是渔场形成的决定性因素，毛里塔尼亚海域为几内亚湾暖流和加纳利寒流冷暖水团的交汇处，携带大量的营养盐，海域初级生产力较高，为浮游生物的繁殖创造了有利条件，因此叶绿素a浓度是形成渔场的另一个重要的海洋环境要素。需要指出的是，虽然关系模型显示叶绿素a浓度是影响短线竹䇲鱼CPUE的一个因素，但由于关系模型中其系数为-0.027 3，因此叶绿素a浓度对短线竹䇲鱼CPUE的影响相对海表面温度影响较小，这是因为叶绿素a虽然对于渔场的形成不可或缺，但形成渔场的关键还要取决于浮游生物的量和规模［23］，而除了叶绿素a浓度外，营养盐、海流等其他海洋要素同样影响着浮游生物的量和规模。牛明香等［5］利用广义可加模型研究时空和环境因子对东南太平洋智利竹䇲鱼渔场的影响时得出相似的结论，认为叶绿素a浓度虽然与资源密度有关联，但对整个模型的贡献程度相对较低，可作为辅助指标。RUMOLO等［10］通过同位素和胃含物分析西西里海峡大西洋竹䇲鱼的空间分布与摄食行为的关系时指出，大西洋竹䇲鱼的空间分布与温度、光合作用的强度有关，此外与盐度也有关系。本文的关系模型中，海表面盐度对短线竹䇲鱼CPUE的影响没有被表征出来，可能是因为调查期间，作业地点的经纬度变化并不大，海表面盐度值上下浮动很小，盐度梯度相较温度梯度低。

本文通过分位数回归建模，得出关系式后建立IHI模型，并结合Arcmap作图（图2，图3）分析，得出毛里塔尼亚短线竹䇲鱼CPUE较高的区域为：17°20′～17°45′W、20°15′～20°50′N；CPUE较低的区域为16°15′～17°40′W、16°50′～18°50′N。这在资源丰度的时空分布上反映为9—12月短线竹䇲鱼逐渐由南向北转移，ARKHIPOV［12］在研究中指出，毛里塔尼亚海域的短线竹䇲鱼等鱼类在寒冷季节会向北聚集，在布朗角附近产卵，距离南部有一定距离，这从侧面证明了本文研究结果的准确性。张鹏等［24］利用FAO报告和相关资料研究中东大西洋中上层小型鱼类资源及其渔业现状时总结得出短线竹䇲鱼最密集分布区在21°00′N附近，这与本文的研究结果在纬度上大体一致，但也存在差异，原因在于研究所使用的样本数据在经纬度上的分布和分辨率不同。

3.3 建议

本文认为IHI模型能较准确地预测短线竹䇲鱼的CPUE。短线竹䇲鱼IHI分布较高和资源丰度较高的区域在17°20′～17°45′W、20°15′～20°50′N，因此建议我国灯光围网船在毛里塔尼亚海域9—12月作业时，尽量使下网区域集中在上述区域，以增加产量，提高经济效益。短线竹䇲鱼不仅是毛里塔尼亚海域小型中上层渔业的重要组成部分，也是金枪鱼类等的主要食物，因此对其养护与管理十分重要。IHI分布较低的区域是16°15′～17°40′W、16°50′～18°50′N，由于短线竹䇲鱼主要摄食浮游生物，因此在9—12月，可以投放少量流木为浮游生物的繁殖创造条件。

3.4 不足与展望

毛里塔尼亚海域的竹䇲鱼主要包括大西洋竹䇲鱼和短线竹䇲鱼，它们有相似的生物学习性，但具体的栖息水层等存在差异，因此在研究竹䇲鱼的时空分布时，深度也是应当考虑的变量之一。本文缺乏短线竹䇲鱼栖息深度的数据，因此空间上未能说明短线竹䇲鱼的垂直分布。此外，西北非海域上升流显著［25］，鱼类的分布还受溶解氧、营养盐的分布等因素影响［26］，因此要更全面地探索短线竹䇲鱼时空分布与海洋环境的关系，还需要考虑上述影响因子。本次调查时间跨度不大，有效数据偏少，因此之后的研究要尽可能增加观测时长和站点数量。另外，如有可能，可以采集各水层的浮游生物量和短线竹䇲鱼的胃含物等，为更深入地研究短线竹䇲鱼的运动模式提供支撑。

致谢：本研究得到了毛里塔尼亚政府的许可，并得到宏东国际（毛里塔尼亚）渔业发展有限公司兰平勇、林兰英、曾志坚、赵松辉、朱国平先生和“福远渔097”号船长及全体船员等大力支持，谨致谢意。