基于空间统计的珠三角地区土地利用生态风险时空特征

王 彬, 梅志雄, 魏俊超, 马 杰

(华南师范大学 地理科学学院, 广州 510631)

当前,我国的城市化进程加快推进,加速了对土地的索取,土地资源利用程度持续增加,对自然和生态环境造成了极大压力,人地矛盾日渐突出,其中以土地利用结构不合理、土地生态风险加剧等问题最为明显[1-3]。研究表明:土地利用是人类活动与自然环境交互过程中的关键部分,与很多生态系统问题密切相关,甚至对区域生态安全起到首要作用[4-5]。不同土地利用方式和过程对区域内生态系统产生累积性的影响,其影响直观地反映在生态系统的结构和组成变化当中,并涉及到许多生态过程,由此产生了生态风险问题[6-7]。在此背景下,为探究区域土地利用变化对生态环境引起的不利影响,寻找降低其影响的方法和途径,开展区域土地利用生态风险分析并提出相应的管控措施,有利于促进土地资源可持续利用和保障生态安全。

土地利用生态风险是指土地利用及其变化造成的生态系统结构破坏及功能降低的可能性[8],近年来受到诸多学者的关注,且日益注重表现其时空分异特性[9-10]。土地利用生态风险在不同尺度上表现出的空间格局差异较大,具有尺度效应[11]。但现有相关研究通常基于一种空间尺度,如以格网作为评价单元(样方尺度)[12-13]或以行政区(省、区县、乡镇)为评价单元[14],或小流域单元尺度等[15]。这些研究一定程度上推动了土地利用生态风险研究的进展,但仍存在如下缺陷:仅限于单一尺度的土地利用生态风险时空变化分析,普遍缺乏多种空间尺度下的土地利用生态风险规律的对比分析而未真正揭示其尺度效应,且主要依赖经验及主观因素选取研究尺度,缺乏客观判断依据[16]。目前土地利用生态风险评价方法主要有基于风险源汇和基于景观格局两类[4,8]。基于风险源汇的方法需要识别风险受体和风险源,涉及变量多且变量间关系复杂,数据获取困难,且对生态风险的空间表达能力不足、未考虑空间异质性;基于景观格局的土地利用生态风险评价方法是以土地利用变化为诱因,直接从景观格局角度定量识别和评估生态本底安全风险[4,17-18]。这类方法所需数据易获取,便于通过空间采样将景观空间结构转化为空间化的生态风险变量,与GIS结合可深入表达生态风险的空间特征、空间异质性及其时序变化[4,7-8,19-20],逐渐成为近年土地利用生态风险评价的主流方法。目前国内外相关研究热点区域主要集中在城市[18]、流域[7,21]、沿海地带[22]以及重点自然保护区[23]等,而针对珠三角地区开展的相关研究较少。仅见叶长盛等[24]评价了1990—2006年珠三角地区土地利用生态风险,但缺乏2006年后的分析,且只基于5 km×5 km格网尺度。

为此,本文再次从土地利用视角,选取城镇化快速发展、土地利用结构剧烈变化导致生态风险较为突出的珠三角地区为研究区,以2000年、2010年、2020年的土地利用数据为基础构建生态风险指数,通过半变异函数工具对比分析样方尺度、行政单元尺度和小流域尺度下土地利用生态风险的空间异质性,并借助空间自相关和GIS空间分析方法揭示其时空变化特征。研究将力图揭示不同尺度下区域土地利用生态风险的空间结构模式、时空规律及其差异表现,探测其尺度效应,弥补现有相关研究的不足,并识别和探寻出需要关注的重点区域和潜在生态风险区,进而根据生态风险等级制定不同土地规划管理措施,为区域土地资源合理利用和规划管理及促进区域生态安全提供决策参考。

1 研究区及数据

珠三角(21°40′—24°30′N,111°20′—115°30′E)位于广东省中南部,毗邻南海,与东南亚隔海相望,总面积约5.48万km2,包括广州、深圳、佛山、惠州、东莞、珠海、中山、肇庆和江门9个地级市。改革开放尤其是新世纪以来,伴随着经济和城镇化的快速发展,该区大量农业用地和生态用地转为建设用地[25]、土地利用结构变化剧烈,严重扰动了生态环境,区域内生态风险问题备受关注[26]。因此,从土地利用视角分析新世纪以来该区土地利用变化所引起的生态风险时空特征,对生态风险防控具有重要理论及实践指导意义。

本文所用的数据主要包括数字高程模型(DEM)和土地利用数据。从地理空间数据云平台(http:∥www.gscloud.cn/)获得2009年30 m分辨率的DEM数据,用于划分分水岭小流域单元。通过全球地理信息公共产品平台(http:∥www.globallandcover.com/)获取2000年、2010年、2020年30 m分辨率的地表覆盖栅格数据,根据研究区实际和本文需要,将该地表覆盖数据重新归类得到100 m分辨率下耕地、林地、草地、建设用地、水体和未利用地6种地类。

2 研究方法

2.1 评价单元划分

选取格网单元、行政区划单元和小流域单元,分别对应样方、行政单元和小流域3种空间尺度,来实现多尺度土地利用生态风险值的空间表达。具体划分情况为:(1) 根据景观生态学原理,通常将研究区内斑块平均面积的2～5倍来确定格网单元的划分尺度[19]。根据研究区土地景观斑块大小实际,选取0.5 km为步长,从0.5 km×0.5 km到5 km×5 km创建10种格网单元,分别得到216 060,53 940,23 998,13 488,8 629,6 025,4 542,3 383,3 061,2 160个格网评价单元;(2) 以研究区下辖569个乡镇为行政区评价单元;(3) 使用ArcGIS的水文分析模块,从DEM提取出3 415个小分水岭作为小流域评价单元。

2.2 土地利用生态风险指数

为综合表达不同景观类型在外界影响下产生的潜在危害性,参考前人研究[6,8,27],选择土地利用类型面积比例、景观干扰度指数和景观脆弱度指数来构建土地利用生态风险指数,计算公式为:

(1)

式中:ERIk为单元k的土地利用生态风险指数;Ak和Aki分别为单元k的面积(km2)和单元k内地类i的面积(km2);Ei和Fi分别为i地类的景观干扰度指数和景观脆弱度指数。

2.2.1 景观干扰度指数 景观干扰度指数Ei表示景观要素受到外部环境干扰的程度[28-29],可表达为:

Ei=aCi+bSi+cUi

(2)

式中:Ci,Si,Ui分别为i地类的景观破碎度指数、分离度指数和优势度指数;a,b,c为权重,参考相关研究[28]并结合专家咨询分别赋值为0.5,0.3,0.2。

(1) 景观破碎度指数(Ci)指某一地类在外界干扰下从均匀的整体状态转变为异质的破碎状态,其值越大,单元内部越不稳定,可由i地类的斑块数与单元面积的比值表示。

(2) 景观分离度指数(Si)指某地类的斑块单元在空间上的分离程度,其值越高,表示该地类的空间散布程度越大,容易受到外界干扰。其计算公式为:

(3)

(4)

(5)

式中:Si为i地类的景观分离度指数;Di为i地类的距离指标;Pi为i地类的面积比值;A和Ai分别为单元面积(km2)和单元内i地类的面积(km2);ni为单元内i地类的斑块数。

(3) 景观优势度指数(Ui)指某一地类对景观整体的影响程度,其值越大,表明该地类对景观格局的影响程度越大,可表达为:

(6)

式中:Ui为i地类的景观优势度;n为单元内总斑块数;ni,Pi含义同上。

2.2.2 景观脆弱度指数 景观脆弱度指数Fi表示景观要素结构受到外部因素扰动时的易损性,一般通过专家咨询法并结合实际对不同地类赋值并归一化获得。参考相关研究[20,28]和专家意见,并进行归一化得到6种地类的脆弱度指数如下:耕地0.190 5,林地0.095 2,草地0.142 9,建设用地0.047 2,水体0.238 1,未利用地0.285 7。

2.3 半变异函数

半变异函数是地统计学中研究区域化变量空间异质性的重要方法,能够表述变量的结构性和随机性变化特征[30]。其计算公式为:

(7)

式中:γ(h)为半变异函数;h为样本点间隔距离即步长;N(h)为步长为h时的样本点对数量;z(xi)和z(xi+h)分别为位置xi和(xi+h)处的土地利用生态风险值。

半变异函数有3个重要拟合参数:(1) 块金值C0,表示区域化变量的不确定性;(2) 基台值C+C0(包括块金值C0和结构方差C),为半变异函数保持平稳时对应的数值,描述变量整体变异情况;(3) 变程A,表示变量保持空间自相关性的距离。常用的半变异函数模型有指数模型、球状模型和高斯模型等,可通过模型拟合后的决定系数R2和残差RSS大小对比来确定最优拟合模型[30]。

2.4 空间自相关方法

空间自相关用来探测空间单元与邻域单元在相同属性值上的关联程度,并揭示其空间异质和空间集聚,以检验区域化变量的空间结构关系,包括全局和局部空间自相关,常用Moran′sI指数表示[31]。

全局Moran′sI反映区域单元观测值在整个研究区的空间自相关特征[32],其表达式为:

(8)

式中:Xi为i单元的观测值;Wij为行标准化的空间权重矩阵;I值范围落在[-1,1],在定义的显著性水平下,当I显著为正时,单元观测值呈空间集聚分布;当I显著为负时,单元观测值呈空间散布状态;当I值为0时,表明单元观测值为独立随机分布。

全局空间自相关虽可探测区域中具有同一属性单元的总体集聚程度,但无法表示局域单元的空间集聚特征。而局部空间自相关不仅可探测空间异质,也可准确反映出集聚区的具体位置和范围[31]。局部Moran′sI表达式如下:

(9)

式中:S为各单元观测值的标准差。

局部Moran′sI计算结果可划分为高—高集聚、低—低集聚、高—低集聚和低—高集聚4种类型。高—高(低—低)型即某单元与邻域单元间表现为正空间自相关,高(低)值在局域空间呈现集聚现象;高—低(低—高)型即某单元与邻域单元间表现为负空间自相关,高(低)值单元被邻域的低(高)值包围,是空间异质区域[31]。

3 结果与分析

3.1 多尺度比较及适宜评价尺度的确定

根据公式(1)—(6),借助Fragstats软件计算得到选择的不同空间尺度下2000年、2010年、2020年各评价单元的土地利用生态风险值,将其作为样本点数据进行后续研究。以2000年土地利用生态风险指数计算结果为例,使用GS+软件计算得到其半变异函数拟合结果(表1)并对比不同尺度的参数变化特征,以选取最适宜研究尺度。

表1 2000年不同尺度下土地利用生态风险指数半变异函数模型拟合结果Table 1 Land use ecological risk index semi-variance function model fitting results at different scales in 2000

由表1可知,0.5 km样方尺度下,土地利用生态风险指数的C0和C+C0均最小,说明该尺度下土地利用生态风险的随机性较小;随着样方尺度评价单元从0.5 km到5 km变化,C0和C+C0逐渐增加,说明土地利用生态风险的随机性逐渐增大,人类活动或自然灾害等外部干扰因素即偶然因素造成的空间异质性逐渐增强。行政单元尺度下的C0和C+C0均最大,说明该尺度下土地利用生态风险的随机性特征比其他尺度都明显。0.5～3 km样方尺度下,其结构方差/基台值即C/(C+C0)呈小幅波动,但均维持在0.82以上。根据地统计学原理,这说明结构性因素引起的空间异质性在样本间总变异中起主导作用,也反映了土地利用生态风险的空间自相关性较强[16]。而3 km～5 km样方尺度下,C/(C+C0)总体呈快速下降态势,说明结构性因素引起的空间异质性逐渐减弱,随机因素的影响程度加深,土地利用生态风险的随机性变异增强。小流域尺度的C/(C+C0)达到0.839,行政单元尺度的C/(C+C0)也较高。从不同尺度拟合结果的决定系数R2来看,0.5～2.5 km样方尺度下的R2总体上比3.5～5 km样方尺度下的R2小,而3 km样方尺度下的R2为0.829,为所有尺度中的最大值;小流域尺度下的R2为0.617,拟合效果尚可;行政单元尺度下的R2为0.453,拟合效果较差。因此,3 km样方尺度是研究区2000年土地利用生态风险分析的最佳尺度。

以同样的方法分别对2010年、2020年的土地利用生态风险进行适宜评价尺度分析,均得到3 km样方尺度为最佳评价尺度。但考虑到小流域尺度的半方差函数拟合效果也较好以及为了保持研究区地理表面特征,最终选择3 km样方尺度和小流域尺度作为本文的分析尺度。

3.2 适宜尺度下土地利用生态风险的空间变异分析

由表2可知,3 km样方尺度和小流域尺度下2000年、2010年、2020年半变异函数最优拟合模型均为指数模型。这两个尺度下各年拟合结果的决定系数R2均较高,且残差RSS均较小,表明各年拟合效果均较好。观察各年最佳拟合参数中的块金值C0,不论是3 km样方尺度还是小流域尺度,均为2010年C0值最小、2020年C0值最大,表明2020年土地利用生态风险的空间异质性受到随机因素的影响程度较强。两个尺度下基台值在2000—2020年均逐渐降低,表明土地利用生态风险的空间变异程度有所减弱,C+C0变化不大说明结构趋于稳定。两个尺度下结构方差和基台值的比值C/(C+C0)均为2010年最大,且2000年和2010年的C/(C+C0)均达80%以上,表明研究区土地利用生态风险的空间异质性受随机部分因素影响较小,以结构性因素影响为主,也反映该区域土地利用生态风险的空间自相关性较强;2010—2020年,两个尺度的C/(C+C0)均明显下降,说明随机因素对土地利用生态风险的影响逐渐增强。两个尺度下2000年变程A均最大,表明该年研究区土地利用生态风险的空间相关性的作用范围较大。

表2 2000-2020年土地利用生态风险指数半变异函数模型拟合结果Table 2 Land use ecological risk index semi-variance function model fitting results for 2000-2020

3.3 土地利用生态风险空间自相关分析

3.3.1 全局空间自相关分析 根据公式(8)计算各年两个尺度下土地利用生态风险的全局Moran′sI值和相关指标(表3)。由表3可知,两个尺度下各年的全局Moran′sI值均显著为正,说明2000—2020年珠三角土地利用生态风险整体呈现正空间集聚模式,即高—高值集聚、低—低值集聚。3 km样方尺度下全局空间正相关性整体比小流域尺度下更强;相较于2000年、2010年两个尺度下的全局Moran′sI值小幅上升,说明2000—2010年两个尺度下的这种正空间集聚程度均略微增强。但2020年两个尺度下的全局Moran′sI值又略有下降,说明2010—2020年两个尺度下的这种正空间集聚程度趋于减弱。

表3 3 km样方尺度和小流域尺度下土地利用生态风险的全局Moran′s I和相关指标Table 3 Global Moran′s I and related indicators of land use ecological risk at the 3 km quadratic scale and small watershed scale

3.3.2 局部空间自相关分析 利用公式(9)计算各年两个尺度下各单元土地利用生态风险值的局部Moran′sI值,然后借助ArcGIS软件空间化得到局部空间集聚格局演化图(图1)。研究期内,两个尺度下研究区大部分评价单元土地利用生态风险的局部空间结构变化较小,高—高集聚和低—低集聚单元相对较多且表现相对突出。高—高集聚区主要集中在较高风险和高风险区域,其主要地类为耕地和水域等,但其面积随着时间推移有减少趋势,说明土地利用生态风险有一定程度的减弱趋向;低—低集聚区主要集中分布在研究区中部,其土地利用类型较单一、多为建设用地,同时该集聚区有向外扩散趋向;低—高集聚区零星分布于高—高集聚区周围,主要地类为耕地和草地等;高—低集聚区零星分布于低—低集聚区周围,主要地类为建设用地和耕地,建设用地的景观脆弱度水平较低、内部结构稳定,对减轻土地利用生态风险值有一定作用。低—高和高—低集聚区零星分布且无多大变化,说明研究期内区域土地利用生态风险的空间异质区域并未发生剧烈的变动,区域内部土地利用导致的生态风险不均衡性虽然总体表现不强烈,但依然存在。

图1 3 km样方尺度和小流域尺度下土地利用生态风险的局部空间集聚格局演化Fig. 1 Evolution of localized spatial agglomeration patterns of land use ecological risks at the 3 km quadratic scale and small watershed scale

3.4 土地利用生态风险的时空变化分析

前文已证明3 km样方尺度和小流域尺度下3个年份的土地利用生态风险指数均具有显著的空间自相关性,因此可利用克里金法进行空间化插值。由于两个尺度下划分的评价单元大小不同,生态风险指数值计算结果相差较大,但两个尺度的计算结果是相对的而不影响生态风险等级对比分析。为方便对比各年土地利用生态风险空间分布状况,采用自然断裂点法将两个尺度下的插值结果均划分为5个等级,其中3 km样方尺度下插值结果划分为低风险区(ERI<0.065)、较低风险区(0.065≤ERI<0.086)、中风险区(0.086≤ERI<0.112)、较高风险区(0.112≤ERI<0.150)、高风险区(0.150≤ERI<0.230),小流域尺度下的等级划分为低风险区(ERI<0.006 2)、较低风险区(0.006 2≤ERI<0.008 3)、中风险区(0.008 3≤ERI<0.010 6)、较高风险区(0.010 6≤ERI<0.014)、高风险区(0.014≤ERI<0.024)。最终得到3个年份土地利用生态风险等级空间格局图(图2)。

图2 2000-2020年3 km样方尺度和小流域尺度下的研究区土地利用生态风险空间格局Fig. 2 Spatial patterns of land use ecological risks in the study area at the 3 km quadratic scale and small watershed scale in 2000-2020

从图2可以看出,3 km样方尺度和小流域尺度下的土地利用生态风险值的空间分布状况比较相似,都随着时间的推移,总体风险程度有所降低。研究区大多都以低风险和较低风险等级为主,主要分布在珠三角的中部、西部和北部地区,包括广州市中心城区、东莞、深圳、肇庆和佛山等地区;低风险区的土地利用类型主要是景观脆弱度较低的建设用地,且有向其周边扩张的趋势;较低风险区主要分布在研究区边缘的山地和丘陵地区,且集中在有大片林地的区域;中风险区散布于整个研究区,其范围从2000—2020年逐渐缩小,以耕地、林地和草地类型为主;较高风险区和高风险区分布范围相对较少,且两个尺度下均有收缩趋势,主要分布在顺德区、南沙区、从化区、博罗县、中山市北部以及珠海市南部的斗门区和金湾区等地。研究区中部的较低和低风险区域范围向外扩张,面积增加,与建设用地集中分布区域有所吻合,而建设用地从碎片化逐渐演变为有序规模模式,能够降低损失程度,对于降低生态风险有较大作用。研究区内较高和高风险区的主要土地利用类型以耕地、草地、水域和碎片化的建设用地等为主,由于耕地和水域的景观脆弱度较大,且耕地、草地与零碎小斑块的建设用地交错分布在一起会造成景观破碎度增加,同时这些区域的人类活动较为频繁,其原有景观的完整性遭到人为破坏,受到的外界干扰程度较大,导致生态风险值相对较高。

由表4可知,3 km样方尺度下2000年和2010年以较低风险区和中风险区为主,面积占比分别为70.67%,69.79%;2020年较低风险区和低风险区占主导,面积占比达到75.88%。2000—2020年,土地利用生态风险值变化总体趋于降低,表现在较低和低风险区范围变大,其中低风险区的面积增加幅度达到了8.91%;而中等风险以上区域范围均呈缩小态势,中风险区面积减少10.06%,较高风险区面积减少4.51%。小流域尺度下,2000年、2010年、2020年研究区主导风险区均为较低风险和中风险区,面积占比均达67%以上,最高达到72.98%。2000—2020年,土地利用生态风险值变化同样趋于降低,最主要表现为低风险区面积增幅达到了13.63%,中风险区面积减少8.08%,较高风险区面积减少4.97%;2000—2010年的土地利用生态风险变化幅度比2010—2020年变化幅度较小,尤其是低风险区面积增加了11.11%,整体表现出风险降低的趋势。研究区整体的土地利用生态风险值有所降低,原因在于珠三角地区加快产业升级转换,以节约集约用地为原则,大力开展“三旧”改造工程,以提高土地利用集约程度,同时落实耕地保护与生态建设有机结合,优化土地利用结构和布局,促进区域经济生态协调发展。

表4 2000-2020年3 km样方尺度和小流域尺度下的土地利用生态风险等级面积和比例Table 4 Area and proportion of land-use ecological risk classes at the 3 km quadratic scale and small watershed scale in 2000-2020

总体上,两个尺度下研究期内土地利用生态风险均逐渐减弱,且生态风险等级在2000—2010年变化幅度较小而2010—2020年变化幅度较大,其中低风险区的面积增加幅度最大,而中风险区面积减少幅度最大。

4 讨论与结论

(1) 研究期内,珠三角地区土地利用生态风险分析的最佳尺度是3 km样方尺度,小流域尺度也较为合适。这两个尺度下,土地利用生态风险的空间异质性均以结构性因素影响为主且存在较强的空间自相关性,而在2020年随机因素引起的空间变异程度略有加强。

(2) 总体上研究区土地利用生态风险具正空间集聚特征,两个尺度下2010年的土地利用生态风险全局正空间集聚程度均比其他两年略强;而局部空间集聚格局变化不大,高—高集聚和低—低集聚单元较多且表现相对突出,分别集中在高风险区和低风险区位置;低—高集聚和高—低集聚区零星分布且变化较小。

(3) 研究期内研究区土地利用生态风险以低风险和较低风险等级为主,整体上逐渐减弱,并且2010—2020年生态风险减少幅度较大。两个尺度下生态风险空间分布格局较为相似,低风险和较低风险区范围均逐渐扩大,以集中分布的建设用地和林地两种地类为主;中风险、较高风险和高风险区范围均逐渐缩小,以耕地、草地和水域为主。由于珠三角地区持续推出了一系列优化土地利用结构政策及生态保护措施,该地生态风险得到一定的控制。

基于以上结论,本文的政策启示表现在:对于低等级土地利用生态风险区,应持续全面推进节约集约建设用地政策,并加强实施重要生态功能保护区和生态安全屏障功能区等生态修复工作;对于中、高等级土地利用生态风险区,应适当调整土地利用管理措施,提高土地集约利用水平,以增强生态系统结构稳定性。同时应合理配置生产、生活和生态空间区域,在加强耕地保护、守好耕地红线的前提下巩固退耕还林等政策成果,并划定生态保护区和生态缓冲区等措施来缓解土地利用生态风险压力,优化国土空间布局以促进区域可持续发展。

本文利用多个景观指数构建土地利用生态风险评价模型,借助空间统计方法和GIS空间分析技术,揭示了珠三角地区土地利用生态风险的时空变化特征,并得到需要关注的重点区域和潜在的风险区域,研究结果可为区域国土空间优化利用和生态风险管控指明方向。但是,本文仍存在一些后续需进一步完善的地方:首先,本文主要从土地利用视角探讨生态风险,而区域生态风险演替过程较为复杂,其影响因素众多,未来可考虑将风险源等相关生态因素纳入研究中,使生态风险研究更丰富和全面;其次,在划分不同土地利用生态风险等级时主要参考前人的方法,目前还没有统一的标准来确定生态风险等级,对于不同区域的研究结果可比性不强,这需要进一步深入研究。