扁舵鲣 (Auxis thazard) 隶属于金枪鱼科、舵鲣属，为暖水性中上层集群洄游的小型金枪鱼类^[1]，广泛分布于热带和亚热带海区，在中国的东海和南海均有分布。南海蕴藏着丰富的金枪鱼资源，最新研究表明，以扁舵鲣为主的南海外海的小型金枪鱼类资源量为85万吨，但捕捞量很少，属于轻度开发的鱼类资源，具有较好的渔业开发前景^[2]。在当前近海渔业资源逐渐枯竭的情况下^[3-5]，对南海扁舵鲣这类小型金枪鱼的开发是转移近海捕捞压力的方法之一。

鱼类资源的时空分布具有明显的空间异质性特征，可采用地统计学方法进行研究。相较于经典统计学方法，地统计学能弥补其随机独立性假设的局限性，在基于空间自相关的计算建模、估值分析方面具有明显优势^[6]。1985年地统计学方法首次应用于渔业领域，用于估计生物量，而后在鱼类资源丰度和空间异质性评估等方面应用较为广泛^[7]。目前国内外学者已对南海扁舵鲣的形态鉴定^[8]、生物学特性^[1]、食性分析^[9]、遗传结构和遗传多样性评价^[10]、资源评估及灯光罩网调查^[11-12]等方面做了较多研究，而有关其时空分布的研究却鲜有报道，目前仅对闽中、闽东渔场扁舵鲣的时空分布与温盐关系进行了探讨。本文利用2016—2017年的南海灯光罩网调查资料，采用地统计学方法进行普通克里金插值，对南海扁舵鲣资源密度指数空间异质性特征进行分析，探索其时空分布变化规律，为合理开发和保护我国南海扁舵鲣资源提供科学依据。

1.   材料与方法

1.1   数据来源

渔业数据来自2016—2017年4个航次的灯光罩网调查数据 (图1)，调查船总吨位421 t，长43.6 m，配备金属卤化物集鱼灯 (1 kW)，所用罩网的主尺寸为281.60 m×81.76 m，网衣最大网目35 mm，网囊最小网目20 mm，沉纲配重2816 kg。调查期间每晚19点开灯300盏，开灯后约3 h进行放网起网作业，并进行渔获物统计，记录内容包括调查时间、经纬度、渔获物种类、体质量等基本信息。4个航次的调查时间依次为2016年4月和10—11月、2017年4—5月和8—9月。

图  1  南海渔业资源调查站点图

Fig. 1  Survey station of fishery resources in South China Sea

下载: 全尺寸图片 幻灯片

对开灯数量和捕捞量进行归一化处理，作为计算单位努力捕捞渔获量 (CPUE) 的基础。处理公式为：

式中：$ {i} $ 为站号；X_i-nor为第i个站点的开灯数量/捕捞量进行归一化；X_i为第i个站点的捕捞量/开灯数量；X_min为该航次的最小开灯数量/捕捞量；X_max为该航次的最大开灯数量/捕捞量。

CPUE作为表征相对资源丰度的指标^[13]。计算公式为：

式中：$ {C}_{i} $ 为归一化后第 $ {i} $ 站的捕捞量；${L_i} $ 为归一化后第 $ {i} $ 站的开灯数量。

1.2   方法

结合经典统计学和地统计学方法分析南海扁舵鲣分布的空间异质性，地统计学以区域化变量理论为基础，要求数据必须满足正态分布，对不符合正态分布的数据可通过对数、平方根、反正弦平方根、倒数和Cox-Box等转换，以达到地统计学分析的要求^[14-15]。该方法由变异函数和克里金插值两部分组成，主要是对特定空间内的变量进行变异性建模，并基于此模型进行变异性估值^[16]。变异函数是地统计学特有的工具，常用的理论模型有球状模型、高斯模型和指数模型。其中，球状模型的一般公式为：

高斯模型的一般公式为：

指数模型的一般公式为：

式中：$ \text{γ}\text{(}\text{h}\text{)} $为变异函数；$ \text{h} $为步长，有效步长为观测步长，与组数的乘积一般不超过研究区域最大距离的一半^[17-18]，本文中的步长设置为样点间的平均距离^[19]，采用弦距离法计算每一个要素的质心到最近的要素的距离之和的平均值；为对变异函数进行生态学解释，可利用理论模型对变异函数进行拟合，拟合效果可用决定系数、残差、标准平均值、均方根、平均标准误差等指标来检测，从而获得最佳拟合模型^[19-20]。

采用Garrison^[21]的分布重心法，对各航次站点的经纬度以CPUE为权重进行加权平均，确定南海扁舵鲣资源密度指数重心变化的轨迹。公式为：

式中：$ \overline{\rm {lon}} $、$ \overline{\rm {lat}} $分别为重心分布的平均经度和纬度；${{\rm {lon}}_{i}} $，$ {{\rm {lat}}_{i}} $分别为第$ i $站位的经度和纬度；$ { {D}}_{i} $ 扁舵鲣在 $ i $ 站位的CPUE；$ {n} $为站位数。不同航次南海扁舵鲣资源重心分析结果的观测误差可能会使得不同随机采样调查下的结果出现偏差。因此采用Bootstrap自助抽样分析分布重心的置信区间，该方法是对样本进行重复抽样1000次，构建分布重心的置信区间，以此来比较南海扁舵鲣不同航次分布中心变化的显著性^[22]。

普通克里金插值法既考虑了空间变量的相关性，又兼顾有效数据的数量，优势明显^[16]。故本文采用考虑了变量空间相关性的普通克里金法进行插值, 对不同年份扁舵鲣CPUE进行插值时，全部选取移除一阶趋势进行处理^[23]。

采用W-S检验进行正态性检验，GS+9.0软件进行变异函数的计算和模型的拟合，调查站点图和克里金插值图由Arcgis 10.2软件绘制，底图数据来源于全国地理信息资源目录服务系统（http://www.webmap.cn/）。

2.   结果

2.1   数据检验与常规统计分析

对各航次南海扁舵鲣资源CPUE数据采用W-S拟合优度检验数据的正态性，发现数据不符合正态分布，进行对数转换后数据均满足正态性要求 (表1)。

表  1  各航次数据W-S正态性检验

Table  1  W-S normality test of each voyage

航次 Voyage	季节 Season	P	转换后P P for conversion
1	春	0.002	0.884
2	秋	0.046	0.521
3	春	0.000	0.091
4	夏	0.000	0.153
注：P>0.05，该组数据具备正态特征。 Note: P>0.05. The set of data has normality.

下载: 导出CSV 

| 显示表格

对4个航次CPUE数据进行常规统计计算，得到数据分布的各项基本特征参数 (表2)，各航次偏度介于1.218 0~3.235 0，均大于0，分布形态均为右偏；除第2航次外，其余航次峰度值均大于3，呈尖峰分布，低资源密度指数的海域所占比重较大，高资源密度指数海域较少；变异系数CV介于0.879 7~2.034 2，除了第2航次为中等变异情况 (0.1<CV<1)，其余均属于强变异程度 (CV>1)，渔场资源密度指数差异较大；另外，从各年际和季节来看，资源密度指数也存在差异。年际差异程度依次为2017年 (CV均值1.814 1)>2016年 (CV均值1.236 5)；季节差异程度依次为春季 (CV均值1.813 8)>夏季  (CV均值1.594 0)>秋季(CV均值0.879 7)。

表  2  扁舵鲣调查数据基本统计参数

Table  2  Basic statistical parameters of survey data of A. thazard

航次 Voyage	最小值 Min.	最大值 Max.	均值 Mean	标准偏差 Standard error	方差 Variance	偏度 Skewness	峰度 Kurtosis	CV=S/m
1	0.003 7	1.500 0	0.326 8	0.520 7	0.271 0	2.074 0	4.153 0	1.593 3
2	0.019 7	1.000 0	0.349 2	0.307 2	0.094 0	1.218 0	0.582 0	0.879 7
3	0.001 4	1.000 0	0.120 7	0.245 6	0.060 0	3.235 0	10.654 0	2.034 2
4	0.000 1	1.000 0	0.183 4	0.292 3	0.085 0	2.060 0	3.547 0	1.594 0

下载: 导出CSV 

| 显示表格

2.2   南海扁舵鲣空间异质性分析

基于南海扁舵鲣空间异质性特征，建立变异函数模型获得相关理论参数 (表3)。结果表明，球状模型在各航次空间异质性结构中发生频率较高；各航次变异函数块金值介于0.000 1~0.035 0，差异较小，各航次CPUE空间异质性受随机因子的影响程度较弱；基台值介于0.250 2~1.230 0，变化趋势与块金值相同，均为先下降后增加；一般采用块金系数度量样本空间相关性，比值表示随机因子引起的空间异质性中自相关部分占系统总变异的比例，各航次块金系数介于0.000 1~0.077 0，属于较强的空间自相关性，在南海扁舵鲣空间分布中，结构性成分起主要作用；变程与资源密度指数和分布范围有关，反映资源密度指数大小，密度较高的集群变程一般较小，各航次变程值介于1.073 9~2.410 0，无明显变化，差异较小。

表  3  各航次扁舵鲣资源变异函数参数

Table  3  Variation function parameters of A. thazard resources in each voyage

航次 Voyage	最优模型 Optimum model	块金值 Nugget	基台值 Sill	变程 Range	块金系数 Nugget/Sill
1	球状模型	0.013 0	0.810 0	1.870 0	0.016 0
2	高斯模型	0.000 1	0.250 2	1.073 9	<0.000 1
3	球状模型	0.035 0	0.452 0	2.410 0	0.077 0
4	球状模型	0.006 0	1.230 0	2.090 0	0.005 0

下载: 导出CSV 

| 显示表格

2.3   南海扁舵鲣资源密度指数的时空变化

基于各航次扁舵鲣空间异质性结构分析扁舵鲣洄游的分布特征 (图2)，各航次渔场分布存在一定差别，但均呈现出明显的片状或斑块状特征，且其洄游分布路线大致呈西南—东北走向。第1航次有两个明显的资源密度指数高值区，主要分布在10°N—12°N海域；第2航次的资源密度指数高值区在10°N—14°N的西部海域，由西向东递减，有两个明显的资源密度指数低值区；第3、第4航次的资源密度指数高值区均位于东北部海域，由东北向西南递减，其中第4航次在西南部海域形成了一个次高值区，变化梯度较第3航次更加和缓。

图  2  各航次南海扁舵鲣空间异质性结构分布图

注： a—d依次对应 1—4 航次。

Fig. 2  Distribution of spatial heterogeneity of A. thazard in each voyage

Note: a−d correspond to Voyages 1−4.

下载: 全尺寸图片 幻灯片

各航次南海扁舵鲣CPUE重心均分布在10°N—12°N，多靠近岛礁附近。以第1航次扁舵鲣CPUE的分布中心为起点 (图3)，第2航次向西北移动至12°N附近，第3航次略向西南移动，但仍偏北，第4航次又略向西北移动，变化幅度较小。同时，南海扁舵鲣CPUE重心航次间分布差异不显著 (P>0.05)，相对集中在调查海域的中部偏东南 (表4)。

图  3  各航次南海扁舵鲣CPUE重心移动轨迹

Fig. 3  Migration trajectory of center of gravity of CPUE of A. thazard

下载: 全尺寸图片 幻灯片

表  4  各航次南海扁舵鲣CPUE重心的置信区间 (95％)

Table  4  Confidence interval for center of gravity of A. thazard of each voyage (95%)

航次 Voyage	纬度 Longitude/(°N)	置信区间 Confidence interval/(°N)	经度 Latitude/(°E)	置信区间 Confidence interval/(°E)
1	10.565 9	[9.978 4, 11.629 7]	115.362 5	[112.624 7, 115.279 2]
2	11.702 4	[10.711 8, 12.283 1]	113.272 6	[113.008 2, 115.207 4]
3	11.495 8	[11.140 2, 12.613 5]	114.671 4	[112.921 6, 114.777 5]
4	11.936 8	[11.095 3, 12.400 5]	114.144 3	[112.552 8, 113.921 2]

下载: 导出CSV 

| 显示表格

3.   讨论

3.1   空间异质性特征

本研究采用地统计方法分析了各航次南海扁舵鲣空间分布的异质性结构，发现其布局总体以低密度海域为主，并具有聚集性特征。空间变异函数的种类揭示了鱼群依赖于当前环境的聚集分布程度^[24-25]，南海扁舵鲣空间异质性格局以球状函数为主，表现出结构性和聚集性较好的空间格局，集聚程度较为明显。从变程来看，第4航次变程范围较大，表明其空间自相关尺度较广，其余航次的变程范围均小于4，说明其分布格局呈聚集分布，这与扁舵鲣为中上层集群洄游的鱼类属性有关，从空间变异函数类型和变程范围两方面均验证了扁舵鲣的聚集分布。此外，从CPUE均值和基台值的关系来看，在一定条件下，资源密度指数较高的航次其空间异质性也高，该结论也在许多已有研究中得到了验证^[26-27]。

3.2   渔场时空格局

南海扁舵鲣的资源丰度及其分布受到时间、空间和海洋环境等多种因素的影响。从4个航次看，各季节资源密度指数差异较大，依次为夏季>春季>秋季，这与一般统计分析结论相同^[28]。3—7月为南海扁舵鲣的繁殖期，种群规模增长迅速，故夏季资源密度指数最大。此外，扁舵鲣还具有明显的西南-东北洄游特征。结合其生物学特性^[1,29]发现，4—6月为南海扁舵鲣的产卵期，洄游至北部海域产卵繁殖，后进行索饵洄游至西南部海域^[30]。从海洋环境看，夏季西南季风控制南海，季风吹动表层海水，使得表层洋流由西南向东流动，表层丰富的营养物质和浮游生物随洋流一起迁移，洋流向大陆架和大陆坡输送丰富的营养物质和浮游生物^[31]，扁舵鲣向西南方向索饵洄游，资源密度指数较高。此外，水层深度也会影响扁舵鲣分布，从图2知，扁舵鲣多栖息于近岸浅水域^[30]，渔场重心靠近岛屿。从海面到海底都有较为充分的阳光透射，使浮游生物在光合作用下迅速地进行繁殖，给扁舵鲣饵料以丰富的营养物质，形成高生产力海区^[32-33]。且浅水水域可在风浪、潮汐、对流的影响下，水体混合充分，底层补充到上层，物质循环快，初级生产力高，整个水体营养好^[34]。扁舵鲣多捕获于近岸浅水海域，不仅与扁舵鲣自身的栖息分布有关，也与渔业生产作业有着直接的关系^[33]。

从各航次南海扁舵鲣CPUE分布重心95％的经纬度置信区间 (Bootstrap法，表4) 可知，不同航次的CPUE重心经纬度置信区间多有重叠，且大部分航次重心均位于置信区间内。

除自身的生物学特性外，海洋环境也会影响鱼类的分布和洄游。由于扁舵鲣具有较强的游泳能力，有能力进行长距离迁移^[28]，所以可明显地观察到其分布重心在厄尔尼诺-南方涛动现象 (El Niño-Southern Oscillation, ENSO)、西太平洋暖池、太平洋十年际气候振动等海洋-大气相互作用系统影响下的变动情况^[35]。通过对比各时期的SSTA3.4指数 (数据来自美国NOAA气候预报中心网站，http://www.cpc.noaa.gov，图4)，发现第1航次处于强厄尔尼诺期，第2航次处于弱拉尼娜期，第3、第4航次均处于正常时期。当ENSO事件发生时，直接引起热带太平洋海域水温的大规模变化，对具有洄游能力的扁舵鲣来说，ENSO引起的水温结构变化则会影响其空间分布和洄游。南海扁舵鲣渔场重心在厄尔尼诺期偏东，拉尼娜期偏西，正常时期位于中间 (图3)，这与Lehodey等^[36]对中西太平洋鲣鱼 (Katsuwonus pelamis) 的研究结果相似，当厄尔尼诺现象发生时，鲣鱼群体整体向东迁移约4 000 km；当拉尼娜现象发生时，则反向迁移4 000 km。但扁舵鲣的迁移范围没有鲣鱼那么大，推测是因为鲣鱼具有高度洄游能力，扁舵鲣的洄游能力不如鲣鱼强。此外，推测ENSO事件的强度与持续时间也会影响扁舵鲣群体的迁移范围。

图  4  各航次Nino 3.4 分析

Fig. 4  Nino 3.4 index analysis of each voyage

下载: 全尺寸图片 幻灯片

3.3   插值方法选择

空间插值是基于已知数据点对研究区域进行预测的方法，在渔业资源领域应用广泛^[37-38]。通常可被分为确定性插值和不确定性插值 (地统计插值)，确定性插值 (如反距离权重法、自然邻域法、趋势面法和样条函数法等) 是以样本点和插值点之间的距离为权重进行加权平均，地统计方法 (如克里金插值法) 是以包含空间自相关的统计学方法为基础，不仅可预测表面，还可对预测结果的准确性提供某种度量。

不同插值方法对前提假设和原始数据的要求不同，常用的反距离权重法适用于均匀分布的点^[39]，其优点是直观高效，发生各向异性时，会考虑方向权重，但研究表明反距离权重法只考虑插值点与样本点之间的距离，依赖于前人的经验，易受观测点数据集的影响^[40]；普通克里金插值法要求数据呈正态分布，其优点为从变量自身特点出发，考虑观测点的整体空间分布情况，可对插值误差做出理论估计，且能给出估计精度，具有平衡性，结果更精确，更符合实际^[37]，但区域化变量离散性太强时，估值不够精确，且只限于单变量在空间分布的特征研究^[18]。本研究采用的普通克里金插值方法要求区域变量满足二阶平稳或本征假设，既考虑了平稳范围的大小，又兼顾有效数据量，是一种折中方案，与其他插值方法相比，普通克里金插值法适用范围广、计算简单，更加符合实际情况。此外，在进行空间插值前，需要对渔业数据进行检查，降低异常情况的影响，其方法有很多，如划分小渔区等，通常越精细的空间尺度越能体现原始数据的空间特点^[41]；还可以利用3倍平均值±标准差、箱图等方法寻找和剔除异常值^[42]。

3.4   展望

本研究由于原始调查站点密度较小，从一些周围的点到极值点过程中ln(CPUE) 渐进不显著，导致内插结果受极值点的影响明显。今后，在调查站点设计环节应充分考虑野外实际情况；考虑积累更长时间序列数据，适当扩大采样范围和采样密度，在扁舵鲣栖息范围内多设置站点，使采样数据更具代表性。