基于地统计与空间自相关的山区草地资源空间分布格局及变化特征分析

孙晓莉，袁 磊，赵 然，和万荣，田淑静，郭朝元

(1.云南国土资源职业学院测绘地理信息学院,云南 昆明 652501; 2.云南师范大学信息学院, 云南 昆明 650500；3. 徐州工程学院管理学院，江苏 徐州 221000)

【研究意义】草地资源是陆地生态系统的重要组成部分，在水土流失防治、生态环境保护等方面起着举足轻重的作用[1-2]。草地资源不仅是发展国民经济的物质基础，还具有维护生态环境的天然屏障功能[3]。近年来，受经济发展、人口增加等因素影响，中国草地资源退化严重[4]，草地资源面积减少。但随着退耕还草、退林还草等草地保护政策的实施，部分区域草地资源面积出现了增加趋势，针对草地资源的这一变化情况，国内外学者开展了大量研究工作。【前人研究进展】主要集中在草地资源类型退化时空变化特征及其原因分析[5-8]、微生物对草地资源数量、质量的影响[9-10]、基于遥感数据的生物量反演[11]、草地植被覆盖变化[12]等方面，传统的研究侧重于草地资源类型、草地资源质量等方面的较多。【本研究切入点】基于地统计和空间自相关模型从县级、乡(镇)级、村级区域单元3种不同尺度定量地对山区草地资源空间分布格局及变化特征研究较少。【拟解决的关键问题】以云南省会泽县为研究区域，从县级、乡(镇)级、村级区域单元不同尺度，基于OpenGeoda、ArcGIS等软件，从草地资源动态变化度、草地资源集中度指数、核密度分析、空间自相关模型对2009-2016年草地资源空间分布格局及变化特征进行分析，为草地资源清查、生态环境保护以及全面深化草原生态文明体制改革提供数据支撑和决策依据。

1 材料与方法

1.1 研究区概况

会泽县位于云南省东北部、金沙江东岸、曲靖市西北部，东经103°03′～103°55′、北纬25°48′～27°04′，地处滇东高原与黔西高原结合部，位于乌蒙山系主峰地段，地势南高北低，呈阶梯状下降，境内最高海拔大海梁子牯牛寨4017.3 m，最低海拔小江四级电站695 m，相对高差3322.3 m。境内草地资源丰富，主要包括山地草甸类、灌木丛类等。会泽县共有21个乡(镇)，2009年土地总面积为588 912.70 hm2，其中草地资源75 928.38 hm2，占研究区土地总面积的12.89 %。

1.2 数据来源与处理

数据来源：2009、2016年草地资源数据来源于遥感影像，图像空间分辨率为0.5 m，主要采用目视解译、监督分类和外业核查的方式获取草地资源图斑；行政区划等数据来源于会泽县第二次全国土地调查数据库。文中主要基于ArcGIS10.0和OpenGeoda软件对数据进行处理为下文研究提供数据基础。

图1 研究区位置示意图Fig.1 The location of research area

1.3 研究方法

1.3.1 草地资源动态变化度 草地资源动态变化度用于定量描述研究区域草地资源面积数量和空间上的变化情况以及在一定的监测时间段内变化的剧烈程度。草地资源动态变化度计算模型如下[13]：

式中：K为草地资源动态变化度，Sb为某区域监测结束时期的草地资源面积，Sa为某区域监测开始时期草地资源面积，T为监测时间段。

1.3.2 草地资源集中度指数 草地资源集中度指数的概念来源于地理集中度指数[14]，是指草地资源在研究区域空间地域上集中程度的测算指标，主要用来分析草地资源在地域空间上的分布情况。草地资源集中度指数计算模型如下[15-16]：

式中：G为草地资源集中度指数，xi为空间区域草地资源分布面积，T为研究区草地资源总面积，n为研究区空间区域个数。

1.3.3 核密度分析 核密度分析属于一种非参数密度估计的统计方法，其是热点和冷区识别与分析的一种有益的探测性方法[17]。在草地资源空间分布格局及变化特征分析中，核密度分析主要是通过输入研究区域的草地资源数据来测算整个研究区域草地资源图斑的集聚情况。假定x1,x2,x3.....xn为独立同分布f的n个样本，设其概率密度函数为f，则核密度分析的计算模型如下[18]：

式中：fn为草地资源地块核密度计算值，n为观测草地资源地块的数量，k为核密度函数，x～xi为测算草地资源地块x到样本草地资源地块xi处的距离，h为核密度计算带宽或平滑参数。

1.3.4 空间自相关 空间自相关是指一些变量在同一个分布区内的观测数据之间潜在的相互依赖性。而Moran’sI是衡量空间自相关集聚程度的基本测度，其包括全局空间自相关和局部空间自相关，在进行空间自相关分析时应首先确定空间权重矩阵。

(1)全局空间自相关：全局空间自相关(Global Moran’sI)是全面测度研究区域地理要素属性值聚集或离散的程度。Global Moran’sI计算模型如下[14]：

标准化Z值经常用于检查全局空间自相关指数的显著性水平，其计算模型如下：

式中：E(I)为自相关性的期望值，VAR(I)为方差，Zscore表示标准化统计量的阈值。

主要选取Global Moran’sI来测算各村级区域单元草地资源面积比例(即草地资源面积占所在村级区域单元面积的比例)的全局空间自相关程度。

(2)局部空间自相关：局部空间自相关(Local Moran’sI)指数可以更好的分析研究区域内各村级区域单元在空间位置上与各自邻近单元属性的相关性，其可以使用图形的形式更直观地表示研究区域各村级区域单元的空间集聚情况，并根据LISA聚集图分析研究区域村级区域单元的聚集或离散的空间位置分布，LISA图的结果常分为HH(高-高型)、LL(低-低型)、HL(高-低型)、LH(低-高型)以及NN(非显著相关型)，其中HH和LL型均表现为空间正相关类型，HL和LH型均表现为空间负相关型，NN则表示无明显的聚集或离散特征。局部空间自相关计算模型如下[19-20]：

主要选取局部 Moran’sI来测算各村级区域单元草地资源面积比例(即草地资源面积占所在村级区域单元面积的比例)的局部空间自相关程度。若Moran’sI>0,则说明草地分布存在空间集聚，差异较小；反之，若Moran’sI<0,则说明草地分布存在空间显著性差异。

2 结果与分析

2.1 草地资源区域分布及变化特征

2.1.1 草地资源面积变化分析 会泽县草地资源

图2 研究区草地资源集中度指数雷达图Fig.2 Radar map of the research area grassland resources concentration index

动态变化度在2009-2016年全县草地资源数量呈现增加趋势，草地资源动态变化度为3.98 %，共增加21 134.51 hm2。研究区域各个乡(镇)草地资源动态变化度最高是田坝乡为20.01 %、其次是雨碌乡(14.24 %)、新街乡(13.41 %)、纸厂乡(12.09 %)、大井镇(11.31 %)。全县21个乡(镇)中，2009-2016年间草地资源增加的乡(镇)有16个，减少的有5个，减少的乡(镇)中动态变化度较大的是马路乡(-2.93 %)、矿山镇(-2.18 %)和老厂乡(-1.85 %)。各个乡(镇)草地资源面积的变化有一定的差异。

2.1.2 草地资源空间分布变化特征 2009年草地资源分布较多的乡(镇)为中西部的大海乡、待补镇、金钟镇，北部的迤车镇和火红乡，南部的田坝乡，而中部的者海镇、大井镇、新街乡则分布的较少，整体呈现出“南北两端高、中间低”的空间分布特征。2016年草地资源分布较多的乡(镇)为大海乡、待补镇、田坝乡、金钟镇等，而位于研究区域中上方的老厂乡、大桥乡、乐业镇、矿山镇则分布较少，与2009年相比，整体呈现“南北两端高、中间低”的空间分布特征更为明显。

如图2所示，全县草地资源集中度2009年为7.70，2016年为7.98，2009-2016年增加了0.28。2016与2009年相比大多数乡(镇)草地资源集中度指数变动幅度较小，其中增加较明显的乡(镇)主要有田坝乡、老厂乡、矿山镇和娜姑镇等；减少较为明显的乡(镇)主要有驾车乡、者海镇、马路乡和五星乡等。从2009和2016年草地资源集中度指数空间分布来看，高值区主要集中在北部的纸厂乡、西北部的五星乡和大桥乡、中部的新街乡、东南部的鲁纳乡以及东北部的矿山镇区域，整体上呈现出“西北、中东部区域高，其他区域稍低”的空间分布特征。

2.2 草地资源核密度空间分布及变化特征

借助ArcGIS10.0软件，首先将2009和2016年的草地资源图斑通过“要素转点”工具，将草地资源图斑(面数据)转换为点数据，然后使用“核密度分析”功能进行草地资源核密度分析，并将结果采用自然间距分类法划分为不同的区间，如图3所示。2009年草地资源地块核密度计算结果，共划分了0～2、2～4、4～6.52块/km23个不同的等级区间。2016年草地资源地块核密度计算结果，共分为了0～2、2～4、4～5.11块/km23个不同的等级区间。2009、2016年研究区草地资源空间分布核密度测算最大值分别为6.52和5.11块/km2。从核密度计算结果来看，部分区域在2009-2016年期间草地资源地块空间分布数量呈现减少的趋势。2009年核密度计算结果高值区域集中在大桥乡的中西部、五星乡的东部、雨碌乡的西南部以及鲁纳乡的西部，2016年高值区域则集中在五星乡的东部。2009-2016年核密度计算高值区域发生明显的减少；研究区域草地资源核密度分析中值区域(2～4块/km2)：2009年主要分布在除研究区两侧区域外的其他区域，主要集中在西北部的大桥乡、五星乡、乐业镇，中部区域的新街乡、雨碌乡和鲁纳乡，南部区域的田坝乡等区域；2016年主要集中在五星乡的竹箐村、铅厂村、大坪子村、红石岩村、野猪冲村、干松林村、石龙村，少量分布在纸厂乡的纸厂村、灯草塘村和浑水塘村，鲁纳乡的窝坡村和哈克村，田坝乡的白岩村、田坝村等区域。中值区域，除五星乡和田坝乡变化不大外，大桥乡、乐业镇、雨碌乡、新街乡、上村乡等发生了不同程度的减少；低值区域(0～2块/km2)，2009-2016年出现了由外围向中心扩张的趋势，但整体变化不太明显。

图3 研究区草地资源核密度计算空间分布Fig.3 Spatial distribution of grassland resources kernel density of reseatch area

2016与2009年相比草地资源地块分布密度较为分散，除五星乡外，无明显集中区域，且研究区两侧草地资源核密度明显低于其他区域，整体呈现出以五星乡为中心的西北部区域较集中，其他区域较分散的情况。从草地资源图斑分布情况来看，2009-2016年期间草地资源面积增加一部分来源于退耕还草和退林还草，空间分布相对较为分散，且增加的草地资源图斑约45 %分布于坡度25°以上的区域；草地资源面积增加仅仅是单个图斑面积增加，草地图斑数量却在减少，导致2009-2016年草地资源核密度测算值出现了减少的情况。

2.3 草地资源分布空间关联性分析

2.3.1 全局空间自相关 借助ArcGIS10.0软件，首先计算出各个村级区域单元的草地资源面积比例，然后再基于相同字段使用空间连接进行关联，分别得出2009和2016年草地资源面积比例图层数据；然后再使用OpenGeoda软件对研究区2009和2016年各村级区域单元草地资源面积比例值进行全局Moran’sI指数，并对其进行显著性检验，结果如表1所示。2009和2016年研究区各村级区域单元草地资源面积比例值全局Moran’sI指数均>0且对比标准化Zscore(2.58)与显著水平的临界值，2009和2016年的Zscore分别为13.63、12.88，均明显大于检验临界值，并通过了显著性水平a=0.05的检验，充分说明研究区域各村级区域单元的草地资源面积比例值呈现出显著的空间自相关性，表现出显著的集聚分布情况，即研究区各村级区域单元草地资源面积比例值高的村级区域单元其周围区域草地资源面积比例值也高，反之，各村级区域单元草地资源面积比例值低的村级区域单元其周围区域草地资源面积比例值也低。

表1 研究区2009和2016年各村级区域单元草地资源面积比例值全局Moran’s I指数Table 1 Global Moran's I index of proportion of grassland resources area at village level in 2009 and 2016

图4 2009和2016年各村级区域单元草地资源面积比例值局部Moran’s I散点图Fig.4 Scatter-plots of local Moran’s I for the proportion of grassland resources area at village level in 2009 and 2016

2.3.2 局部空间自相关 研究区各村级区域单元草地资源面积比例值的全局空间自相关分析结果显示2009和2016年草地资源虽然在空间分布上呈现出集聚状态，但却出现了减弱的趋势。从草地资源动态变化度、草地资源集中度指数、草地资源核密度等空间分布情况来看，各个乡(镇)内部存在较大的差异性。因此，通过局部空间自相关分析来研究区域内516个村级区域单元草地资源面积比例的空间自相关水平的差异性显得尤为重要。

首先利用OpenGeoda软件获取研究区2009和2016年各村级区域单元草地资源面积比例的局部Moran’sI指数，然后以其为横轴，各村级区域单元草地资源面积比例的局部Moran’sI指数相对应的空间滞后向量值为纵轴，绘制出研究区2009和2016年各村级区域单元草地资源面积比例观测值的局部Moran’sI指数散点图，如图4所示。

从图4可以得出，如果横坐标>0，则说明空间单元为草地资源面积比例观测值较高的村级区域单元，反之则为较低的村级区域单元；如果纵坐标>0，则空间单元的邻近单元为草地资源面积比例观测值较高的村级区域单元，反之则为较低的村级区域单元。根据图4的局部Moran’sI散点图结果，将其划分为四个象限。其中第I象限为正相关的“高-高”类型区、第Ⅱ象限为负相关的“低-高”类型区、第Ⅲ象限为正相关的“低-低”类型区、第Ⅳ象限为负相关的“高-低”类型区。

从研究区2009-2016年间草地资源面积比例局部Moran’sI指数散点图中可以得出：①2009和2016年位于第I象限的“高-高”类型区、第Ⅲ象限“低-低”类型区的村级区域单元占比较重，其中2009和2016年的个数分别为370、357个，即位于这2个象限的草地资源面积比例较高(或较低)的区域呈现出较明显的空间集聚分布特征，且位于“低-低”类型区的村级区域单元个数高于“高-高”类型区的个数；2009和2016年位于第Ⅱ象限的“低-高”类型区、第Ⅳ象限“高-低”类型区的村级区域单元占比较轻，其中2009和2016年的个数分别为146、159个，且位于“低-高”类型区的村级区域单元个数高于“高-低”类型区的个数。②2016与2009年相比，位于各个类型区的村级区域单元发生了不同程度的变化，其中“高-高”、“低-低”类型区的个数有所减少，说明呈现空间集聚分布的空间单元在减少，呈现空间分散特征的空间单元在增加。2016与2009年相比位于第I象限的“高-高”类型区单元减少了7个，但是从空间分布情况来看，除大海乡附近空间集聚分布情况未发生较大变化外，其他区域的空间分布发生了较大的变化，由研究区域的东北部转变为了中东部和南部区域；2016与2009年相比位于第Ⅱ象限的“低-高”类型区单元增加了8个，从空间分布情况来看除增加了研究区的南部区域外，其余区域未发生较大的变化；2016与2009年相比位于第Ⅲ象限“低-低”类型区单元减少了6个，空间分布未发生较大的变化；2016与2009年相比位于第Ⅳ象限“高-低”类型区增加了6个，空间分布增加了北部区域，其他区域未发生较大变化。③从整体来看，研究区域2016与2009年相比，尽管位于第I和Ⅲ象限的村级区域单元个数在减少，整体上呈现集聚分布的趋势在减弱，但总体变化不太大。

为了更好、更加直观的表示研究区域各村级区域单元草地资源面积比例自相关水平差异程度以及变化趋势等，在a=0.05的显著性水平下，使用OpenGeoda软件，并借助于ArcGIS10.0软件获得研究区各村级区域单元草地资源面积比例LISA集聚图，如图5所示。以期更好的表示空间单元与相邻单元同质和异质的局部特征变化及其地理分布情况。

从LISA集聚图中可以看出“高-高”、“低-低”、“高-低”、“低-高”类型区所表现的特征主要有如下。

(1)“高-高”类型区：该类型区是指其自身属于草地资源面积比例高的村级区域单元，周围相邻单元也是草地资源面积比例高的村级区域单元，属于局部均质分布的高值集聚区。2009和2016年“高-高”类型区的村级区域单元数量分别为52、40个。2009年主要分布于大海乡，娜姑镇的东南部，金钟镇的西部，待补镇的中东部、中西部，驾车乡的西北部，上村乡的中部，迤车镇中东部，马路乡的西南部，火红乡东部，矿山镇北部等区域；2016年主要分布于大海乡、娜姑镇东南部、待补镇西南部、驾车乡西北部、大井镇的中西部和西南部、者海镇南部、雨碌乡中部、田坝乡东南部等区域。2016与2009年相比，大海乡周边区域仍保持“高-高”集聚区不变，主要是由于大海草山旅游景区一定程度上促进了草地资源保护，该区域地势平缓，草地集中连片，且草地资源十分丰富。而其他区域呈现出由研究区东北部(迤车镇、矿山镇等)向中东部(大井镇、雨碌乡等)以及南部(田坝乡)偏移的趋势。另外，受到退草还林以及退草还耕等因素的影响，也使得2009-2016年期间出现“高-高”集聚区的变化中既有增加也有减少(缩减)情况的出现。

(2)“低-低”类型区：该类型区是指其自身属于草地资源面积比例低的村级区域单元，周围相邻单元也是草地资源面积比例低的村级区域单元，属于局部性空间差异较小的低值集聚区。2009和2016年“低-低”类型区的村级区域单元数量分别为44、42个。2009和2016年“低-低”类型区主要分布于者海镇、新街乡和金钟镇3个乡(镇)的交界处、乐业镇和矿山镇的交界处以及大桥乡的西北部等区域。2016与2009年相比，虽然“低-低”类型区总体位置分布差异不大，但较为明显的出现了缩减和扩大，主要表现为者海镇、新街乡和金钟镇3个乡(镇)交界处出现了缩减情况，新街乡不再有“低-低”类型区，者海镇也出现了一定的缩减；扩大情况出现在以乐业镇东侧为边界，向乐业镇东北部以及矿山镇的西部进行扩大，其他区域变化不大。研究区“低-低”类型区出现缩减和扩大的原因主要是：影像判读时，当种植的树木较小时，影像判读为草地；随着树木的生长，茂密的林地遮挡了草地，便判读为了林地，这样便出现了部分区域草地资源面积的增加和减少，从而影响了研究区“低-低”类型区的缩减和扩大情况。

(3)“高-低”类型区：该类型区是指其自身属于草地资源面积比例高的村级区域单元，周围相邻单元是草地资源面积比例低的村级区域单元，呈现空间负相关性，空间异质性较强。2009和2016年“高-低”类型区的村级区域单元数量分别为2、0个。该类型区分布较少，2009年主要分布在大井镇的井田村和金钟镇的仙龙社区，2016年研究区没有该类型区的分布。该区域出现缩减的原因主要是出现了村庄及建制镇的占用情况导致该区域草地资源面积减少，而邻近单元草地资源面积出现了增加，使得彼此草地资源面积比例对比不再明显。

(4)“低-高”类型区：该类型区是指其自身属于草地资源面积比例低的村级区域单元，周围相邻单元是草地资源面积比例高的村级区域单元，呈现空间负相关性，空间异质性较强。2009和2016年“低-高”类型区的村级区域单元数量分别为10、16个。“低-高”类型区在空间上分布较为分散，2016与2009年相比，“低-高”类型区出现了不同程度的缩减和扩大变化。

3 讨 论

许多学者以县域为研究尺度，从不同草地资源类型的植被特征来研究草地资源空间分布特征[6]、草地资源退化时空分布特征[21]等，对空间分布关联性研究较少。但近年来随着GIS技术的快速发展[22]，运用GIS和相关地统计模型对草地资源空间分布格局进行研究成为一种可能，从不同尺度[县、乡(镇)、村级区域单元]、不同角度(宏观和微观)对研究区草地资源动态变化进行定量分析，突出了草地资源空间分布的动态变化及其空间分布的关联性，不仅对开展草地资源空间分布及变化特征研究具有重要的借鉴意义，而且对传统研究中缺乏草地资源动态变化及其空间分布关联性进行了补充。

本文从数量和空间上对草地资源空间分布格局及变化特征进行研究，但草地资源动态变化是一个长期的过程，空间分布特征及其关联性往往受研究区温度、降雨量、地形条件、耕作半径等因素影响。为了更加精细化、定量化、动态化的表达草地资源的动态变化特征，下一步将基于遥感影像建立实时的草地资源动态变化监测体系[23]，进一步探讨影响草地资源退化/恢复的驱动力因子，不断促进草地资源的可持续利用。

图5 各村级区域单元草地资源面积比例值LISA集聚图Fig.5 LISA plots for the proportion of grassland resources area at village level

4 结 论

从县级、乡(镇)级、村级区域单元3种不同的尺度，借助于地统计与空间自相关模型对研究区2009-2016年草地资源在宏观、中观和微观不同层次的空间集聚和差异性特征进行定量分析，揭示了研究区在2009-2016年期间草地资源数量和空间上的动态变化。结果表明草地资源动态变化度较大的为田坝乡，草地资源相对较为稳定的区域为大海乡、娜姑镇及待补镇，与当地通过草原生态奖补机制项目的贯彻落实加快该区域草原生态恢复、防止草原退化以及与多部门联合执法杜绝破坏草原现象发生等举措具有密切关系。此外，研究结果可作为当地草地资源保护区划分以及为今后生态建设和可持续发展等提供重要依据。