任何事物或现象和空间位置具有依赖关系，空间位置越近，则事物或现象越相近，空间位置越远，则越相异。海洋鱼类通常以个体、种群、群落的形式分布在特定的时间和空间上，因此具有高度的时空异质性。摸清鱼类的空间分布格局及其动态对认识人类活动和气候变化背景下鱼类种群动力学及对渔业管理政策的制定具有重要意义。

在事物空间特征研究方面，探索性空间数据分析(exploratory spatial data analysis，ESDA)方法得到越来越广泛的应用。ESDA以空间关联性测度为核心，其目的为描述与显示对象的空间分布，发现异常观测值，揭示空间联系、簇聚以及其它异质性的空间模式^[1-2]，从庞杂的数据中发现它们在空间上的联系及变化规律。作为ESDA技术的核心，空间自相关分析自1950年Moran等^[3]提出以来，国外学者较早地将其应用到各个研究领域。Jumars^[4]利用自相关分析方法，对半深海底栖息生物群落的纵向和横向空间结构进行了研究。Benzie^[5]则对广泛分布的海洋种类刺海星(Acanthaster planci)的空间自相关进行了分析。Ye等^[6]利用ESDA方法对1960~1995年芝加哥的自杀率和社会经济因子之间的关系进行了研究。国内学者利用ESDA分析方法进行了一系列研究，如区域经济发展演变及其差异^[7]、城镇群体空间结构^[8]以及土壤元素空间分布^[9]等。但在渔业上的应用少有报道，冯永玖等^[10]基于ESDA方法对西北太平洋柔鱼资源空间热点区域及其变动情况进行了研究。对于黄海鳀分布年际变动特征的研究，目前只对其总体空间分布进行了分析^[11]，而从不同站位资源密度之间的相关性探究其内部空间格局的研究尚未见报道。

鳀(Engraulis japonicus)作为黄、东海生态系统的关键种，自20世纪90年代大规模开发以来其产量显著下降^[12]，空间分布也发生了显著变化^[11]，但其空间格局特征及变化尚不明晰。本文利用ESDA的空间自相关分析方法并结合GIS技术，分析黄海中南部鳀越冬群体的空间分布格局及其动态，并探讨其空间格局与海表温度的关系，以期为气候变化和人类活动影响下的黄海鳀资源管理提供科学依据。

1   材料与方法

1.1   数据来源

鳀资源密度数据来自“北斗”号调查船2002~2015年期间在黄海中南部(121°E~125.5°E、32.5°N~ 37°N)的底拖网调查，取样时间均为1月中上旬，站位布设以经、纬度0.5°×0.5°为一个渔区，各年间站位设置一致。调查按预设站位昼夜进行，各站计划拖网1 h左右，平均拖速3 n mile/h，实际拖网时间根据调查时的海况确定。当场记录渔获重量和尾数，各站位的拖网时间均标准化为1 h，以每小时获得的渔获重量代表鳀相对资源密度。本研究的目的是分析鳀资源的空间分布特征，为去除零渔获站位空间关联对分析结果的影响，分析数据只包含有鳀鱼捕获的采样站位。

对应时段的水温数据源自美国航空航天局(NASA)反演的全球卫星遥感三级数据产品，时间分辨率为1个月，空间分辨率为4 km(http://oceancolor.gsfc.nasa.gov)，利用IDL读取原始数据提取研究区的水温，取1个渔区范围内的水温平均值作为该渔区水温。

1.2   数据预处理

变量满足正态分布是进行空间自相关分析的前提，对每个年份的鳀相对资源密度进行检验，发现原始数据不服从正态分布，易导致空间自相关分析结果存在偏差，必须进行数据变换。本研究对原始数据取自然对数，并进行Shapiro-Wilk tests检验，所有年份的Sig.值均大于0.05，满足数据正态分布要求。因此，文中进行空间自相关分析的鳀相对资源密度数据均为取自然对数后的数据。

1.3   空间自相关性分析

1.3.1   全局空间自相关

空间自相关是检验某一要素的属性值与其相邻空间点上的属性值是否显著相关联的重要指标，统计内容包括空间变量的空间位置和属性^[13]。全局空间自相关用于整体定量描述研究变量的空间分布模式，探索某一属性的数据在区域中的整体分布格局，常用的测度指标是Moran’s I指数，其计算公式^[14]如下：

式中:n是参与分析的鳀捕获站位数；x_i、x_j分别为站位i和站位j的鳀相对资源密度；x是所有鳀捕获站位相对资源密度的平均值；w_ij是n*n的空间权重矩阵值，用来定义i和j之间的连接。本研究采用K=4的最近邻矩阵构建空间权重矩阵，为每个站位选择距离最近的4个站位为邻域。

Moran’s I取值范围为[-1, 1]。在给定显著性水平时，I值大于零，代表空间集聚，表示某站位与其“邻居”之间的相似性大于差异性；小于零，代表空间孤立；等于零，表示不存在空间自相关性，绝对值越大表明变量的空间结构性越强。

1.3.2   局部空间自相关

局部空间自相关揭示局部直至每个空间单元的空间自相关性质，通常用LISA(local indicators of spatial association)指数表示，由Anselin^[15]于1995年提出，其本质上是将Moran’s I分解到各个区域单元，对于某个空间单元i，计算公式为：

各变量的含义同式(1)。LISA指数用于评估局部空间模式的显著性，识别热点区域，进一步研究其空间分布规律。

局部空间自相关分析结果的另一种表达形式为Moran散点图。散点图的横轴为变量x的所有观测值，纵轴为空间滞后向量(w_x)的所有取值，每个区域观测值的空间滞后通过标准化的空间权重矩阵进行定义，其值为该区域周围邻居观测值的加权平均。散点图将空间格局分为四种类型，即高高(H-H)、低高(L-H)、低低(L-L)和高低(H-L)类型，其中H-H和L-L代表空间集聚，而L-H和H-L代表空间孤立。H-H位于第一象限，表示鳀高资源密度站位被其它高资源密度站位所包围；L-L位于第三象限，表示鳀低资源密度站位被其它低资源密度站位所包围；L-H位于第二象限，表示鳀低资源密度站位被其它高资源密度站位所包围；H-L位于第四象限，表示鳀高资源密度站位被其它低资源密度站位所包围。此处H和L是相对于区域总体的平均水平而言，不同年份没有可比性。

1.4   空间分布格局与水温的关系

利用鳀空间格局分布图和海表温度图叠加，直观、定性分析黄海中南部越冬鳀空间分布格局和海表温度的关系；根据越冬鳀适宜海表温度范围及其代表性海表温度值，利用ArcGIS栅格计算器计算121°E~126°E、32.5°N~37°N范围内代表性等温线不同年份的平均纬度，利用产量重心法^[16]计算鳀高高类型集聚区资源密度重心纬度，并进行相关性分析，以探求高高类型集聚区在纬度上的变化和海表温度的关系。产量重心法计算公式如下：

式中:X、Y分别为产量重心的经纬度；X_i为渔区i中心点的经度，Y_i为渔区i中心点的纬度；C_i为渔区i的产量；n为渔区总个数。

1.5   分析软件

利用SPSS软件对鳀相对资源密度进行正态分布检验；用Geoda 0.9.5-i软件进行空间自相关性分析；用ArcGIS 9.3进行相关图件的绘制，且由于原始数据为WGS84经纬度坐标，对渔获数据分布图进行投影和坐标转换(Beijing 1954 3 Degree GK CM 123E)，以满足统计分析需要在平面坐标系下进行的要求。

2   结果与讨论

2.1   鳀空间分布总体特征

表 1为2002~2015年黄海中南部鳀越冬群体资源密度取对数后的全局Moran’s I估计值及其显著性。从表中可以看出，整个研究期间，除2002年、2004年外，全局Moran’s I值均为正且大于0.1，表明大多数年份鳀资源密度水平相似(高高或低低)的站位在空间上呈集聚分布，但各年份的空间集聚程度不同，2015年空间集聚性最强，其次为2010年，都存在显著集聚现象。2005、2006年的空间集聚程度相同，其显著性均为10%，呈弱集聚格局。2002和2004年，Moran’s I值接近于零，鳀在空间上随机分布。2002~2015年，黄海中南部越冬鳀空间自相关程度的变化并无年间变动规律。这是因为影响鱼类数量和分布变动的因素非常复杂，主要因子有捕捞、气象、地理和生物等^[17]，多种因子的综合作用导致了鱼类的空间分布特征，而捕捞、气象和生物等因子存在不规律的年间变动，因此越冬鳀的空间自相关程度也呈现不规律的年间变动。

表  1  黄海中南部越冬鳀分布全局自相关分析结果

Tab.  1  Results of global Moran's I for wintering anchovy distribution in central and southern Yellow Sea

下载: 导出CSV 

| 显示表格

2.2   鳀空间分布局部差异

Moran散点图可以用来研究局部的空间不稳定性，进一步区分单元与邻居之间的空间联系形式。图 1是2002~2015年黄海中南部越冬鳀空间分布的Moran散点图。由图可知，2003、2010、2011和2015年，鳀分布站位多数位于第一、三象限内，分别达到66.67%、89.29%、66.67%和76.74%，为正的空间联系；2005、2006年，鳀分布站位位于二、四象限的数量增多，即呈现空间孤立的站位数量增多，但一、三象限的数量高于二、四象限的数量，分布在一、三象限的站位数量分别为62.5%和63.89%，仍呈现空间集聚特征；2002、2004年，位于一、三象限的鳀站位数量和二、四象限的数量基本相等，鳀空间上呈随机分布。

图  1  2002~2015年黄海中南部越冬鳀相对资源密度Moran’s I散点图

Fig.  1  Scatter plot of Moran's I for the relative stock density of wintering anchovy in central and southern Yellow Sea during 2002 and 2015

下载: 全尺寸图片 幻灯片

为了更加清晰了解黄海中南部越冬鳀各种空间关联模式的空间位置，进一步研究其空间分布规律，将各个年份的Moran散点图绘制到研究区域的地图上形成Moran散点地图(图 2)。从图中可以看出，2002和2004年，仅有少数H-H站位聚集分布，L-L站位分散分布，且H-L和L-H站位比例较大。2003年，L-L集聚区主要分布在121.5°E~123°E、34°N~36°N范围内，H-H聚集区主要分布在123.5°E~124.5°E、33.5°N~36°N之间的区域。2005年，L-L聚集区靠沿岸分布，而H-H聚集区主要分布在东部深水区。2006年，H-H聚集区主要分布在122°E~124°E、34°N~35°N以及124°E~124.5°E、35.5°N~36.5°N的区域内，L-L聚集区则在122.5°E~124.5°E、32.5°N~33.5°N范围内分布。2010年，H-H和L-L呈正空间关系的类型都集中成片分布，H-H类型区集中分布在调查区域的南部，为122.5°E~124.5°E、33.5°N~35.5°N范围内，而L-L聚集区主要分布在调查区域的北部，为122.5°E~124.5°E、36°N~37°N区域内。2011年，主要以H-H类型聚集区为主，分布在调查区域的东部，123°E~124°E、33°N~35.5°N范围内。2015年，H-H和L-L类型聚集区成片分布，L-L类型聚集区分布在调查区域的北部近岸水域，H-H类型聚集区分布在东部水域。总体来看，各个年份的空间分布格局变化较大，但高高类型聚集区均分布在调查区的中东部水域，而低低类型聚集区则主要分布在近岸浅水区及北部水域。赵宪勇^[18]的研究表明，黄海鳀的冬季分布有北部和西部浅水、低温水域个体小，中东部深水区个体大的特点。小个体分布区通常对应了低的资源量，陈云龙^[19]的研究也证实了这一点。本论文对越冬鳀空间格局的分析与上述结论基本一致。

图  2  黄海中南部越冬鳀分布格局(Moran散点地图)与等温线分布图

Fig.  2  Distribution pattern of wintering anchovy (Moran scatter map) and SST isotherm in central and southern Yellow Sea

下载: 全尺寸图片 幻灯片

Moran散点地图(图 2)对越冬鳀四种空间类型进行了可视化，但缺乏相应的显著性水平指标。为了进一步探求越冬鳀的空间特征，对2003年、2010年和2015年三个鳀分布显著集聚的年份进行LISA聚类分析，并将聚类结果展布在研究区地图上，获得鳀分布的Moran显著性水平地图(通过p=0.1检验)(图 3)。从图 3可以看出，在相同的显著性水平下，2010和2015年的空间分布格局明显，而2003年的空间分布格局较弱。Moran显著性水平图同时反映了鳀分布的热点区域。结果表明，2003年，鳀低值分布的热点区域明显，近岸浅水区域分布，高值分布的热点区域不明显；2010年，鳀高值分布的热点区域明显，分布在研究区南部，而低值分布的热点区域则主要在北部；2015年，鳀高值和低值分布的热点区域都较为明显，低值热点区域分布在北部近岸水域，而高值热点区域则分布在东部水域。

图  3  黄海中南部越冬鳀相对资源密度的Moran显著性地图

Fig.  3  Moran significance map of the relative stock density of wintering anchovy in central and southern Yellow Sea

下载: 全尺寸图片 幻灯片

2.3   鳀空间分布格局与海表温度的关系

图 2表征了鳀空间分布格局与海表温度的叠合，图中等温线由高值向低值凸出为暖水，由低值向高值凸出为冷水，图中的暖流为黄海暖流。由图可知，高高类型集聚区主要分布在黄海暖流暖水舌及以东区域，空间分布集聚明显的年份(2003、2005、2006、2010、2011和2015)，主要分布在122°E~124.5°E、33.5°N~36.5°N之间，随黄海暖流的年间变动而出现空间分布上的变化。低低类型集聚区的分布在2003、2005和2011年主要在近岸浅水区，2006年则分布在33.5°N以南的区域，2010年分布在36°N以北区域，而2015年分布在近岸浅水区。总体来看，空间分布集聚年份的低低类型区主要分布在冷暖水交汇的区域。水温作为海洋生物生活环境中最基本最重要的因素之一，不仅对鱼类的繁殖、生长和新陈代谢有影响，而且对鱼类的洄游影响也较大^[20]，水温变化影响鱼类的空间分布和集群^[21-22]。前人研究表明，鳀小个体补充群体的适温范围较大个体要低^[23]，与高高类型集聚区分布在黄海暖流暖水舌及以东区域，而低低类型集聚区主要分布在冷暖水交汇的区域这一空间分布格局相一致。

高高类型聚集是黄海越冬鳀空间分布格局的重要类型，能代表鳀资源密度的主要分布趋势。因此，本研究对鳀高高类型分布与水温的关系进行了定量分析，进一步明晰鳀分布格局与水温的关系。由图 4可知，高高类型聚集区主要出现在海表温度为7~12℃的范围内，因此以该范围内最为集中的8℃、9℃、10℃、11℃表层水温等值线代表鳀空间分布的适宜温度因子，分析其平均纬度变动与高高类型集聚区鳀重量重心纬度变化的关系，以研究水温变化对其分布的影响。

图  4  鳀高高类型集聚渔区出现频率与海表温度的关系

Fig.  4  Relationship between fishing area frequency of high-high spatial aggregation areas of anchovy and sea surface temperature

下载: 全尺寸图片 幻灯片

Pearson相关性分析结果如表 2所示，从表中可以看出，分析数据不包含2004年时，所有等温线的平均纬度均与重量重心纬度呈正相关，除9℃等温线外，两者的相关性均显著(p < 0.05)；当分析数据包含2004年时，8℃和9℃等温线平均纬度与重量重心相关性较差，而10℃和11℃等温线平均纬度与重量重心显著相关，这是因为2004年鳀分布为随机分布状态，高高类型无集聚性。

表  2  鳀鱼高高类型集聚区重量重心纬度与等温线平均纬度相关性分析

Tab.  2  Correlation analysis between the latitude of anchovy weight center of high-high spatial aggregation areas and the mean isotherm latitude

下载: 导出CSV 

| 显示表格

从图 5可以看出，10℃和11℃等温线平均纬度在所有时间序列上和高高类型集聚区中鳀重量重心的纬度变化趋势相一致，除2004年外，8℃和9℃等温线平均纬度和重量重心纬度的变化趋势基本一致，该结果和表 2的分析结果相一致，表明高高类型集聚区沿纬度方向上的移动受海水温度变化的影响。黄海中南部越冬鳀空间分布变动与水温的关系研究^[11]表明，黄海越冬鳀空间分布纬度方向上的变动由水温变化决定，本论文对代表鳀主要分布趋势的高高类型空间分布与水温关系的研究结果与其吻合。冬季，黄海暖流对越冬鳀的分布具有重要影响^[24]，黄海暖流变强的年份，暖水舌向北突出，鳀鱼的分布重心北移，反之，分布重心南移。

图  5  等温线平均纬度与鳀鱼高高类型集聚区重量重心纬度的关系

Fig.  5  Relationship between mean latitude of isotherm and the latitude of anchovy weight center of high-high spatial aggregation areas

下载: 全尺寸图片 幻灯片

鱼类的空间分布受诸多因素影响，温度是最基本也是最主要的环境因子。本研究表明，鳀资源密度的高高类型集聚区的空间分布与海表温度存在相关关系，低低类型聚集区的分布也有一定规律性，但对于不同年份，其空间分布自相关程度的差异原因仅利用海表温度尚无法解释。这是因为环境条件对鱼类行为的影响首先取决于鱼类本身的状况，包括鱼类个体的大小，不同生活阶段和生理状况等^[25]，加之捕捞、生物、地理环境等因素的影响^[17]，以上因子在不同年份都有不同变化，同时因子之间的相互作用十分复杂，本论文尚无法进行详细讨论，这也是有待进一步研究的内容。

3   结论

(1) 不同年份，黄海越冬鳀的空间自相关程度不同。在研究的8个年份中，2015年鳀分布存在中等强度的空间自相关，呈明显空间集聚分布，其次为2010年、2003年和2011年；2005年和2006年鳀分布存在正的空间自相关，但相关程度较小；2002年和2004年Moran’s I值约为零，鳀空间分布处于随机状态。

(2) 黄海中南部越冬鳀空间分布格局年间差异较大，但资源密度的高高类型集聚区主要分布在中东部水域122°E~124.5°E、33.5°N~36.5°N范围内，低低类型集聚区则主要在西部近岸浅水区和北部水域分布。

(3) 越冬鳀高高类型集聚区分布在黄海暖流暖水舌及以东区域，而低低类型集聚区主要分布在冷暖水交汇的区域；且不同年份高高类型集聚区纬度方向上的移动受海水温度变化的影响。

(4) 本文利用ESDA和GIS相结合的方法对2002年~2015年黄海中南部越冬鳀的空间分布格局及其与水温的关系进行了研究，获得了很好的结果，表明ESDA方法可以探索渔业资源空间分布特征，为渔业资源空间格局研究提供了新途径。