矿区土地利用生态风险的时空异质性

常小燕,李新举,*,李西灿,郭 鹏,高 峰

1 山东农业大学信息科学与工程学院,泰安 271018 2 山东农业大学资源与环境学院,泰安 271018 3 山东省济宁市国土资源局,济宁 272001

生态风险是指一个种群、生态系统或整个景观的正常功能受外界胁迫,系统内部某些要素或其本身的健康、生产力、遗传结构、经济价值和美学价值发生退化的可能性[1]。1992年美国环境保护署(U.S. Environmental Protection Agency)将生态风险评价(Ecological Risk Assessment,ERA)定义为“评估暴露于一种或多种压力因子后,可能出现或正在出现的负面生态效应的可能性过程”[2]。区域生态风险评价作为生态风险评价的一个分支,强调在区域尺度上描述和评估环境污染、人为活动或自然灾害对生态系统及其组分产生不利作用的可能性及其大小的过程[3]。空间异质性是指某种生态学变量在空间分布上的不均匀性及复杂程度,即系统特征在时间和空间上的复杂性和变异性[4]。

我国区域生态风险评价起步较晚,主要是对一些生态敏感脆弱区[1,5-19]及城市[20-28],以“风险源-风险受体-暴露危害分析-生态终点”这一概念模型,构建综合生态风险评价指标体系,并借助遥感(RS)和地理信息系统(GIS)技术完成研究区生态风险评价[1,3,5-16,20]；也有一些学者通过对遥感影像解译获取土地利用类型数据,运用景观格局指数构建土地利用生态风险指数,借助空间统计学、地统计学等研究方法,评价土地利用变化所带来的区域生态风险[17-19,21-29]。但不足之处主要在于基于概念模型进行土地利用生态风险评价时,对土地利用生态风险的空间异质性缺乏深层次的分析,仅限于生态风险的状态分析；运用地统计学进行土地利用生态风险空间异质性分析时,大多数研究者人为地根据研究区大小及经验选取研究尺度,缺乏科学依据。

从景观生态学角度将土地利用格局与生态风险相结合以评价区域生态环境,结合空间统计学理论进行土地利用生态风险空间关联性及时空变异性研究,有利于揭示区域土地利用生态风险的时空分布格局及空间变异规律,分析生态风险与土地利用类型、人类活动之间的耦合关系,为土地资源的优化配置及开发整治提供理论基础,以促进土地资源的可持续利用。微山县作为山东省的“南大门”,境内有丰富的矿产资源,同时山东省最大的淡水湖—南四湖也坐落其中；本文从微山县自身生态风险的特殊性入手,在现有矿区生态风险评价研究成果的基础上,选取适宜的景观格局指数构建景观意义上的土地利用生态风险指数,结合空间自相关分析及半变异函数理论,选取适宜的研究尺度,进行研究区土地利用生态风险半变异函数模型拟合,在此基础上进行空间插值分析,以探究2000—2016年微山县适宜尺度下土地利用生态风险的时空变化情况及空间异质性,实现矿区土地利用生态风险空间异质性的定量化和可视化表达,以期为后续的矿区生态恢复与重建工作提供理论依据。

1 研究区概况

微山县位于山东省济宁市南部,地处116°34′—117°24′E,34°27′—35°20′N,南北相距120 km,东西相距8—30 km,总面积1779.8 km2；其中湖面面积1266 km2,占全县总面积的三分之二,由北到南依次为南阳湖、独山湖、昭阳湖、微山湖,统称南四湖,占全省淡水量的45%。区内辖10个镇,2个乡和1个县经济开发区(2014年行政区划)。在选取研究范围时,主要选择煤矿较集中、地表塌陷较严重的地区,同时考虑研究区的连贯性,最终选取微山县的11个乡镇作为研究区,总面积1176.86 km2。研究区地理位置及范围见图1。

图1 研究区地理位置及范围Fig.1 The location and scope of the research area

2 数据来源与研究方法

2.1 数据来源与处理

选取2000年8月21日、2005年9月4日、2010年9月18日的Landsat5 TM遥感影像和2016年9月2日的Landsat8 OLI遥感影像作为数据源,其云量均小于2%,空间分辨率为30 m；影像的时间间隔大致为5年左右,以保证矿区土地利用生态风险评估时间粒度的一致性；同时选用降雨量较充沛、植被覆盖度较高的9月份左右的遥感影像,便于矿区土地利用信息的提取。其他辅助数据包括：微山县土地利用现状图、微山县矿区分布图、微山县2014—2016年国民经济和社会发展统计公报。

参照《土地利用现状分类和编码》(GB-T21010—2007),考虑研究区土地的用途、利用方式、覆盖特征和遥感影像的分辨率,将研究区景观类型分为耕地、其他农用地、城乡建设用地、水域、塌陷积水区、滩涂沼泽6类。利用ENVI 5.3,采用监督分类中的支持向量机分类、人工目视判读与决策树分类、NDWI、NDVI相结合的分层分类,得到2000、2005、2010、2016年4期土地利用类型分布图(图2)。在4期分类图上分别随机选取200个检查点,通过实地调查并结合相近年份的土地利用现状图和天地图网站上的高清卫星图,获取检查点的实际土地利用类型,利用ENVI软件计算混淆矩阵和Kappa系数,经验证Kappa系数均在85%以上,能满足本研究的需要。

图2 2000—2016年研究区土地利用类型分布图Fig.2 Distribution maps of land use type in the research area from 2000 to 2016

2.2 采样方法

运用ArcGIS中的Fishnet(渔网)工具创建0.5 km×0.5 km、1 km×1 km、1.5 km×1.5 km、2 km×2 km的正方形规则格网(样区),其样区个数分别为：4710个、1172个、521个、296个,根据公式计算每一个样区的土地利用生态风险指数值,并将其作为样区中心点的值[18-19,23-29],从而实现生态风险值的空间化表达,以更好地对研究区土地利用生态风险的空间异质性进行研究分析。

2.3 土地利用生态风险指数的构建

2.3.1景观干扰度指数Ei

不同景观类型在保护物种、维护生物多样性、保持生态系统结构和功能完整性、促进景观结构自然演替等方面的作用是有差别的,同时会受到外界环境的干扰,而不同景观类型对外界环境的抗干扰能力也是不同的。参考相关文献[17-18,23-24,27-29],选取景观破碎度指数、景观分离度指数和景观优势度指数,通过这三个指数的叠加构建景观干扰度指数Ei,以反映不同景观类型所代表的生态系统受到干扰(主要是人类开发活动)的程度,其表达式为：

Ei=aCi+bSi+cDOi

(1)景观破碎度指数Ci：景观破碎化是由于自然或人为干扰所导致的景观由单一、均质和连续的整体趋向于复杂、异质和不连续的斑块镶嵌体的过程,景观破碎度是景观异质性的重要组成部分,破碎度值越大,表明景观单元内部稳定性越低。公式为：

Ci=ni/A

式中,Ci为景观类型i的破碎度,ni为景观类型i的斑块数,A为景观(一个样区)的总面积。

(2)景观分离度指数Si：是指某一景观类型中不同元素或斑块个体分布的分离程度,值越大,表明景观在地域分布上越分散,景观分布越复杂[17],受到的干扰程度越大。

Si=Di/Pi

式中,Si为景观类型i的分离度,Di为景观类型i的距离指数,Pi为景观类型i的面积指数,Ai为一个样区中景观类型i的面积,ni、A含义同上。

(3)景观优势度指数DOi：景观优势度表征景观结构中某一类型支配景观的程度,反映了该景观类型对景观格局形成和变化影响的大小[17]。公式为：

DOi=(1+景观类型i的密度+景观类型i的面积指数)/3

其中,密度=(一个样区中景观类型i的斑块数/此样区中总的斑块数),景观类型i的面积指数同上。

根据以上公式计算出Ci、Si、DOi指标后,由于量纲不同,需进行归一化处理。a、b、c为各景观指数的权重,且三者相加为1。权重值的大小反映各景观指数解释景观类型受干扰的程度,根据研究区实际情况,并结合前人研究成果[17-18,23-24,26-29],认为破碎度指数最重要,其次为分离度指数和优势度指数,三个指数分别赋值为0.502、0.301、0.197。

2.3.2景观脆弱度指数Fi

景观脆弱度表示不同生态系统的易损性,与其在景观自然演替过程中所处的阶段有关。一般情况下,生态系统处于初级演替阶段时,食物链结构简单、生物多样性指数小,其较为脆弱。在土地利用生态系统中,人类活动是主要的干扰因素之一,土地利用程度不仅反映土地本身的自然属性,同时也反映人为因素与自然环境因素的综合效应。参考相关文献[11,17-18,23-26]并根据研究区实际情况,本研究区6种景观类型,以塌陷积水区最为脆弱,其次是其他农用地,最稳定的是城乡建设用地,6种景观类型分别赋值为：塌陷积水区为6、其他农用地为5、滩涂沼泽为4、耕地为3、水域较稳定为2、城乡建设用地为1。采用反正切函数归一化方法处理[11,17,23-25]之后,得到各景观类型的脆弱度指数Fi值分别为：0.8949、0.8743、0.8440、0.7952、0.7048、0.5。

2.3.3土地利用生态风险指数ERI

基于上述所建立的景观干扰度指数和景观脆弱度指数,结合各景观类型的面积比重构建土地利用生态风险指数,用于描述一个评价单元内综合生态损失度的相对大小,以便通过采样方法将景观空间结构转化为空间化的生态风险变量。土地利用生态风险指数ERI计算公式如下：

式中,ERIk为第k个采样区内土地利用生态风险指数值,Aki为第k个采样区内土地利用类型i面积,Ak为第k个采样区面积,Ei为土地利用类型i的干扰度指数,Fi为相应的脆弱度指数。

2.4 空间统计分析

2.4.1空间自相关分析

空间统计学是以具有地理空间信息特性的事物或现象的空间相互作用及变化规律为研究对象,以区域化变量为基础,研究既具有随机性又具有结构性,或空间相关性和依赖性的一门科学[30-31],其核心是空间自相关性、空间依赖性和空间异质性。

所谓空间自相关性是指空间上越靠近的事物或现象越相似[4]。空间自相关的基本度量是空间自相关系数,可用全局和局部两种指标来度量。

全局空间自相关(Global Spatial Autocorrelation)用于描述区域单元某种属性值的整体分布状况,判断该属性值在空间上是否存在聚集性的特点。常用的指数是全局Moran′sI指数。

全局Moran′sI指数用来评估整个研究区域内所有空间对象是集聚分布、离散分布还是随机分布。其计算公式为[30-35]：

2.4.2地统计分析

(1)半变异函数

本研究借助地统计学中的半变异函数深入探讨景观结构变化等造成的土地利用生态风险的时空变化情况和空间变异规律。半变异函数的数学公式[4,32-36]为：

式中,r(h)是变异函数,h是两样本配对抽样的分隔距离,即步长,Z(xi)和Z(xi+h)分别是随机变量Z在空间位置xi和xi+h上的取值,N(h)为分隔距离为h时的样本点总对数。

以r(h)为纵坐标,h为横坐标做图,得到半方差图(semivariogram)。半方差图中包含3个重要参数：1)块金值(nugget)C0,表示区域化变量在小于抽样尺度时的非连续性变异[35],2)基台值(sill)C+C0,即平稳值,描述变量在研究区域范围中总的空间变异程度,基台值包括块金值(C0)和结构方差(C)两部分。3)变程(range)A,表示极限距离,即某特征在空间上自相关的空间幅度。

(2)克里格插值

如果变异函数和相关分析的结果表明区域化变量存在空间相关性,就可以利用区域化变量的原始数据和变异函数的结构特点,对未采样点的区域化变量的取值进行线性无偏、最优估计[35]。

3 结果与分析

3.1 土地利用类型结构变化特征

由表1及图3可知,17年间研究区耕地面积呈波动式减少,2000—2005年面积减少4604.94 hm2,2005—2010面积又有所增加,至2016年耕地面积为28706.85 hm2；城乡建设用地2000—2005年增幅最大,增加了2830.41 hm2,其后面积呈缓慢增加趋势,2010—2016年间建设用地减少了804.96 hm2,17年间建设用地共增加4589.19 hm2；塌陷积水区2000—2005年面积增幅最大,之后有所缓解,2010—2016年面积减少了132.57 hm2,整个研究期内塌陷积水区面积共增加473.04 hm2；其他农用地和水域面积总体呈增加趋势,滩涂沼泽呈减少趋势。

表1 2000—2016年研究区土地利用类型面积变化

图3 2000—2016年研究区土地利用类型面积变化统计图 Fig.3 Statistic map of land use type area change in the research area from 2000 to 2016

3.2 空间自相关分析

根据空间概念化模型中的“共享边或角”规则建立空间权重矩阵,并进行行标准化,计算四种尺度下土地利用生态风险的全局Moran′sI指数(表2),结合P值和z得分,进行自相关显著性检验[37]。以1 km×1 km尺度下的土地利用生态风险指数值为例,首先进行零假设验证,以正态分布95%置信区间双侧检验阈值1.96为临界值,根据表2,研究期的P值均小于0.05,z得分均大于1.96,说明研究区土地利用生态风险的空间分布不是随机模式,存在空间自相关性。1 km×1 km尺度下2000年、2005年、2010年、2016年的Moran′sI指数分别为0.4576、0.4950、0.4279、0.4677,说明土地利用生态风险的空间分布呈集聚分布模式,具有很强的空间正相关性。

表2 1 km×1 km尺度下全局Moran′s I 指数及验证值

3.3 土地利用生态风险的尺度效应分析

空间异质性依赖于尺度(粒度和幅度),粒度和幅度对空间异质性的测量和理解有重要的影响[4]。为了使建立的半变异函数模型能准确地反映各种尺度上的变化特征,要确定最佳的采样尺度[4,38]。根据土地利用生态风险指数的构成,以2005年土地利用类型分布图为试验数据,选取破碎度指数,并计算研究区Shannon多样性指数值,借助GS+地统计分析软件,进行不同尺度下的破碎度指数和Shannon多样性指数半变异函数模型拟合(表3、表4),以确定土地利用生态风险空间变异的最佳研究尺度。

表3 不同尺度下破碎度指数半变异函数拟合模型参数

根据表3,0.5 km尺度下,破碎度指数的C0为0.0014,C+C0为0.0040,相比较1 km、1.5 km和2 km尺度下的C0和C+C0,其值最小,说明该幅度内尺度效应不明显；从0.5 km到1.5 km,C0从0.0014上升到0.0056,说明人类活动等随机因素引起的空间异质性在逐渐增强[4],1.5 km到2 kmC0从0.0056骤然下降到0.0004,尺度的增加导致破碎度发生了明显变化。C+C0从0.5 km到1 km上升最快,1 km到2 km缓慢上升,尺度的增加会降低空间自相关性。0.5 km到1 kmC/(C+C0)上升最快,之后缓慢上升,也反映了空间相关性在逐渐减弱。综上可知,0.5 km尺度过小,1.5 km到2 km空间自相关性明显减弱,1 km是适宜的研究尺度。

表4 不同尺度下Shannon多样性指数半变异函数拟合模型参数

分析不同尺度下多样性指数的半变异函数模型拟合参数(表4),尺度由0.5 km上升到1.5 km时,C0从0.0228下降到0.0206,到2 km时C0降至最低值0.0183。根据半变异函数的理论,块金值较大,说明较小尺度上的某种过程不可忽视,随机因素引起的空间异质性起主要作用。0.5 km到1 kmC+C0增加至0.0773,1 km到2 kmC+C0由0.0773降低至0.0633,说明随着尺度的增加,多样性指数的变化逐渐减弱。0.5 km到1 kmC/(C+C0)由0.606上升至0.730,说明0.5 km到1 km,多样性指数的空间自相关性逐渐增强；1 km到2 kmC/(C+C0),波动式下降,由结构性因素所引起的空间异质性在不断减弱。因此,1km能较直观的反映多样性指数的空间变异情况。

综上,在进行微山县11个乡镇的土地利用生态风险时空异质性研究时,1 km是最佳研究尺度。

3.4 适宜尺度下土地利用生态风险的半方差分析

在GS+中,分别对4期1 km×1 km的样点数据进行半变异函数值的计算,并进行土地利用生态风险半变异函数模型拟合(表5及图4),其中步长设置为格网的间距1 km,有效滞后距离39 km为采样点最大间距的1/2。

表5 2000年—2016年土地利用生态风险半变异函数模型拟合参数

由表5及图4看出,研究区2000年和2016年土地利用生态风险的最优拟合模型为球状模型,2005年和2010年的最优拟合模型为指数模型。4个时期的土地利用生态风险的空间变异特征主要表现在以下几方面：1)从各个年份最优拟合模型中块金值C0的变化来看,2016年C0最大,表明土地利用生态风险受人类活动等随机因素引起的空间异质性较大,较小尺度上的某种过程不能忽视。2)基台值C+C0,从2000年的0.01000上升至2005年的最大值0.01276,到2016年有所降低,说明土地利用生态风险度空间变异程度随时间的推移在不断增强,其中2005年基台值相对突出,表明这个时段影响土地利用生态风险度的某些因子的空间变异性明显增强,2005年土地利用生态风险总的空间变异程度较高。3)结构方差与基台值的比值C/(C+C0),2000、2005、2010、2016年分别为57.5%、68.8%、67.9%和54.5%,均大于50%,说明气温、降水、地形地貌、土壤类型等空间结构性因素引起的空间异质性占主导地位,从2000年到2005年,其值先逐渐增大至最大值68.8%,说明研究区土地利用生态风险的空间分异从2000年到2005年受结构性因素的影响逐渐加大,生态风险指数在小尺度上的随机变异逐渐被较大尺度上的空间结构性变异所取代。2005—2016年呈缓慢下降趋势,表明随着社会经济的不断发展,人类活动对土地利用干扰程度的不断加深,导致2005年到2016年土地利用生态风险由随机因素引起的空间变异程度在增强。4)变程A表示土地利用生态风险在空间上自相关的空间幅度,在大于该变程的空间范围内,土地利用生态风险的空间自相关性消失,2000年、2005年、2010年、2016年土地利用生态风险的变程分别为13 km、23.16 km、12.75 km、15.19 km,由2005年变程值较大,也可以看出土地利用生态风险值的变化2005年较其他年份剧烈。

图4 2000—2016年研究区土地利用生态风险半变异函数模型拟合曲线图Fig.4 The fitting curve of semivariogram model of land use ecological risk from 2000 to 2016 in the research area

3.5 土地利用生态风险的克里格插值及精度评定

利用GS+完成2000—2016年土地利用生态风险半变异函数最优模型拟合,获取各模型相关参数之后,结合ArcGIS地统计分析模块,在步长1 km下,选择普通克里格插值方法进行空间估值,通过交叉验证进行插值精度评定(表6)。通过验证分析,误差平均值、误差标准平均值最接近于0,均方根预测误差最小,平均标准误差最接近于均方根预测误差,标准均方根预测误差接近于1,预测误差皆符合要求,预测精度较高。

表6 插值预测误差统计表

3.6 土地利用生态风险的时空异质性分析

将插值结果转化为栅格格式,并按照研究区范围进行不规则裁剪。通过ArcGIS数据处理,发现4期土地利用生态风险预测值的范围都介于0.2718—0.5814值域范围内,为了便于比较4期土地利用生态风险的空间分布情况,利用ArcGIS软件中的重分类工具,进行土地利用生态风险等级划分,具体划分5个等级：低生态风险区(ERI<0.32)、中低生态风险区(0.32≤ERI<0.43)、中生态风险区(0.43≤ERI<0.48)、中高生态风险区(0.48≤ERI<0.5)、高生态风险区(ERI≥0.5),最终得到2000、2005、2010和2016年四期土地利用生态风险等级图(图5)。

图5 2000—2016年研究区土地利用生态风险空间分布图Fig.5 Spatial distribution map of land use ecological risk from 2000 to 2016 in the research area

(1)时间变化

从图5及表7、表8可知,2000年研究区以中低生态风险区和中生态风险区为主,占总面积的57.58%,主要分布在微山湖湖区东西两侧,呈集中连片分布,另外独山湖、昭阳湖湖区及周边也以中低生态风险和中生态风险为主,景观类型主要为湖泊水面及其他农用地类型中的台田鱼塘；中高生态风险区和高生态风险区面积相对较少,分别占区域总面积的17.32%和17.36%,呈散列式分布,主要分布于留庄镇东南部、欢城镇北部和中部一小部分、傅村街道中东部、昭阳街道中东部、韩庄镇东部、高楼乡西北部及西南部的带状区域及西平乡的北部,对照同时期的土地利用类型分布图,高生态风险区主要分布在城乡建设用地及工矿企业周边、道路沿线,另外因煤矿塌陷造成零散分布的耕地其生态风险值也较高。低生态风险区面积9125.48 hm2,比例最少,主要分布在微山湖湖心区域。

表7 土地利用生态风险各年面积及比例

2005年以中低生态风险区和中生态风险区为主,占区域总面积的60.66%,其空间分布与2000年的大致一致,面积比例与2000年相比,中生态风险区面积比减少了8.90%,中低生态风险区面积比增加了11.98%。中高生态风险区和高生态风险区面积相比2000年都有所减少,分别减少了8.36%和3.50%,低生态风险区面积19476.69 hm2,占区域总面积的16.51%,相比2000年面积比增加了8.77%。总体上2000—2005年,土地利用生态风险由中、中高、高生态风险向中低、低生态风险转变。

2010年以中低生态风险区和中生态风险区为主,占区域总面积的74.71%,和2005年相比中低生态风险区和中生态风险区面积分别增加了9.06%和4.99%,增幅较小；低生态风险区面积和2005年相比变化不大,基本保持稳定；中高生态风险区和高生态风险区面积相比2005年有所下降,特别是高生态风险区面积降幅较大,下降了10.60%。2005—2010年土地利用生态风险由中高、高生态风险向中、中低生态风险转变。

表8 土地利用生态风险各年面积变化及比例

2016年研究区也以中低生态风险区和中生态风险区为主,占区域总面积的70.51%,和2010年相比,中低生态风险区面积有所减少,减少了14%,但中生态风险区面积增幅较大,增加了9.8%；低生态风险区面积和2010年相比减少了8.51%；中高生态风险区和高生态风险区面积相比2010年又有所升高,特别是中高生态风险区面积增幅较大,增加了9.17%。2010—2016年,土地利用生态风险由低、中低生态风险向中、中高和高生态风险转变。

总体上2000—2016年,土地利用生态风险呈现出由中高、高生态风险向中、中低生态风险转变的特征,其中高生态风险区面积降幅最大,下降了10.56%,中低生态风险区面积增幅最大,增加了7.05%,其次是中生态风险区,面积增加了5.89%,低生态风险区和中高生态风险区面积变化幅度不大。

(2)空间变化

为了更好地分析2000—2016 年间研究区土地利用生态风险的空间异质性演变规律,在ArcGIS中运用Spatial Analyst模块中的Raster Calculator工具,将2000、2005、2010和2016年土地利用生态风险空间分布(图5)进行空间运算,获取2000—2005年、2005—2010年、2010—2016年、2000—2016年的土地利用生态风险空间变化图(图6),将差值运算结果等于0的区域定义为生态风险“稳定”区,等于-1的区域定义为生态风险“改善”区,等于-2、-3、-4的区域定义为生态风险“明显改善”区,等于1的区域定义为生态风险“恶化”区,等于2、3、4的区域定义为生态风险“明显恶化”区。

根据图6及表9,2000—2005年研究区土地利用生态风险程度较低,生态系统较稳定,生态风险稳定区和改善区面积占区域总面积的78.92%,且分布较集中,主要分布在独山湖、昭阳湖、微山湖湖区及周边,景观类型主要为湖泊水域,这部分区域景观连通性较高,人为干扰性较低,生态系统较稳定；生态风险恶化区和明显恶化区分布较分散,主要分布在留庄镇主城区以南新安煤矿周边,欢城镇作为煤矿较集中的乡镇,其生态风险恶化及明显恶化的区域较大,主要集中在蒋庄煤矿、岱庄生建煤矿,柴里煤矿等矿区及周边,另外傅村街道东南部、昭阳街道东北部塌陷积水较严重的区域、赵庙镇孔庄煤矿及周边生态风险恶化程度也较高,生态风险恶化区和明显恶化区主要位于矿区周边及道路沿线,这些区域土地利用人为干扰因素较强且离湖区较近,地面塌陷较严重,景观类型的破碎化程度较高。

图6 2000—2016年研究区土地利用生态风险空间变化分布图Fig.6 Spatial distribution map of land use ecological risk in the research area from 2000 to 2016

2005—2010年生态风险明显恶化区比例降低,稳定区比例较2000—2005年增幅较大,增加了10.52%,且煤矿较集中的欢城镇、傅村街道和昭阳街道生态风险有了大幅改善,主要原因是2005年以来济宁市政府开展了大规模的塌陷地复垦治理工程,微山县作为采煤塌陷较严重的区域,也开展了大规模的复垦治理工程,经过5年多土地复垦,塌陷地的生态环境状况有了一定的改善,耕地经过整治后变得集中连片,塌陷积水区整治后变成规整的鱼塘等水产养殖用地。

表9 土地利用生态风险时空面积变化及比例

2010—2016年生态风险恶化区面积比例大幅增加,主要分布于独山湖、昭阳湖等湖周边的滩涂沼泽,另外欢城镇中东部、傅村街道东部、夏镇街道东北部、韩庄镇东北部等煤矿集中区生态风险恶化区也有成片分布。主要原因是2014年南四湖遭遇了罕见的旱灾,水面和湿地大面积萎缩导致水质缓冲能力差,持续的旱情使独山湖、昭阳湖、微山湖等湖区及周边生态环境遭到严重破坏；同时随着科学技术的发展及采煤机械化程度的提高,矿山的回采率逐步提高,煤炭资源的大规模开采,产生了越来越严重的地面塌陷、矸石堆积、水体污染等生态环境问题。

4 讨论与建议

(1)在进行微山县11个乡镇的土地利用生态风险时空异质性研究时,1 km是最佳研究尺度,在该尺度下,土地利用生态风险的空间分布呈集聚模式,且具有很强的空间正相关性。2000—2005年生态风险指数在小尺度上的随机变异逐渐被较大尺度上的空间结构性变异所取代,2005年土地利用生态风险总的空间变异程度较高,2005—2016年,随着社会经济的不断发展,人类活动对土地利用干扰程度的不断加深,土地利用生态风险由随机因素引起的空间变异程度在增强。

(2)2000—2016年,土地利用生态风险呈现由中高、高生态风险向中、中低生态风险转变的特征。其中高生态风险区面积降幅最大,中低生态风险区面积增幅最大,其次增幅较大的是中生态风险区,低生态风险区和中高生态风险区面积变化幅度不大。17年间土地利用生态风险有了大幅降低,生态风险稳定区和改善区占区域总面积的79.21%,生态风险恶化区和明显恶化区主要分布在东部矿区周边及道路沿线等人类活动较频繁、地面塌陷较严重的区域。

根据以上结论,结合《济宁市采煤塌陷地治理总体规划(2009—2020年)》,在后续的土地利用过程中,宜采取以下治理措施：1)对于东部煤矿集中区,采煤塌陷较严重,这些区域应根据地形、土壤、积水区深浅及地表塌陷的实际状况因地制宜地采取复垦治理措施。对于塌陷较浅的区域,应充分发挥南四湖湖底淤泥养分充足的资源优势,在保护河湖生态平衡的前提下,对沿湖(河)矿井采用引湖填充方法复垦塌陷地,并最大限度的恢复为耕地,增加农用地面积；对于塌陷较深,地面坡度较大的区域,应采取挖深垫浅的治理措施,整理成鱼塘,用于水产养殖业。2)在研究区西南部微山湖沿岸东西两侧,耕地的生态风险恶化程度2010—2016年较高,对此,在保证耕地总量动态平衡的前提下,应积极实行退耕还湖、退耕还湿的治理措施,增加芦苇沼泽湿地的面积,从而保护生物物种的多样性,提高生态系统的稳定程度。

本研究基于景观生态学构建土地利用生态风险指数,但在景观干扰度指数的权重赋值时,没有进一步分析社会、经济因素与景观干扰度指数值间的定量关系,只是概念性的进行综合性度量,所以计算出的土地利用生态风险指数值也是相对的,并不具有绝对性；另外,在进行研究区土地利用生态风险半变异函数模型拟合时,虽然考虑了土地利用生态风险的尺度效应,但缺乏对生态风险的变化方向效应分析,未来可在这两个方面开展进一步的研究。