矿区土地利用生态风险的时空异质性
常小燕1,2 ,李新举2,* ,李西灿1 ,郭 鹏1 ,高 峰3
1 山东农业大学信息科学与工程学院,泰安271018 2 山东农业大学资源与环境学院,泰安271018 3 山东省济宁市国土资源局,济宁272001
摘要: 以微山县煤矿较集中、地面塌陷较严重的11 个乡镇作为研究区,以遥感影像和土地利用现状数据为数据源,获取4 期土地利用类型分布图,根据景观生态学理论构建景观意义上的土地利用生态风险指数,结合空间统计学及地统计学理论,探究2000 —2016 年研究区适宜尺度下土地利用生态风险的时空变化情况及空间异质性。结果表明：1 km×1 km 规则网格是最佳研究尺度,在该尺度下研究区土地利用生态风险的空间分布呈集聚模式,且具有很强的空间正相关性；2000 —2005 年土地利用生态风险指数在小尺度上的随机变异逐渐被较大尺度上的空间结构性变异所取代,2005 年生态风险总的空间变异程度较高,2005 —2016 年随机因素引起的空间变异程度增强；2000 —2016 年土地利用生态风险呈现由中高、高生态风险向中、中低生态风险转变的特征；生态风险恶化区和明显恶化区主要分布在东部矿区周边、道路沿线、地面塌陷较严重的区域,研究区西南部微山湖沿岸东西两侧,耕地的生态风险恶化程度2010 —2016 年较高；17 年间微山湖湖区的生态风险稳定区和改善区大幅增加,生态风险明显下降。
关键词： 土地利用；景观格局；生态风险；尺度；空间自相关；时空异质性
生态风险是指一个种群、生态系统或整个景观的正常功能受外界胁迫,系统内部某些要素或其本身的健康、生产力、遗传结构、经济价值和美学价值发生退化的可能性[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 景观干扰度指数E i
不同景观类型在保护物种、维护生物多样性、保持生态系统结构和功能完整性、促进景观结构自然演替等方面的作用是有差别的,同时会受到外界环境的干扰,而不同景观类型对外界环境的抗干扰能力也是不同的。参考相关文献[17-18,23-24,27-29],选取景观破碎度指数、景观分离度指数和景观优势度指数,通过这三个指数的叠加构建景观干扰度指数E i ,以反映不同景观类型所代表的生态系统受到干扰(主要是人类开发活动)的程度,其表达式为：
E i =aC i +bS i +cDO i
如此复杂而丰富的实践，给鲁迅先生带来的激励、磨练，超乎想象了。鲁迅杂文的识见高、味儿浓，单用“腹有诗书气自华”来解释，一定简单了、片面了、肤浅了——“三杂”，尤其是阅历杂，造就了不可复制的鲁迅先生，成就了鲁迅先生三个伟大！
针对互联网+外贸环境下外贸业务员职业岗位职能的变化，国际贸易实务课程不应仅仅满足传统贸易下职业能力的培养，而应同时满足跨境电子商务交易（B2B、B2C）模式下职业技能的训练。因而，国际贸易实务课程应以培养具有互联网思维和国际贸易知识、同时胜任传统贸易线下交易和跨境电子商务线上交易的销售人才为目标。
C i =n i /A
式中,C i 为景观类型i 的破碎度,n i 为景观类型i 的斑块数,A 为景观(一个样区)的总面积。
(2)景观分离度指数S i ：是指某一景观类型中不同元素或斑块个体分布的分离程度,值越大,表明景观在地域分布上越分散,景观分布越复杂[17],受到的干扰程度越大。
S i =D i /P i
为了使研究结论更为可靠，本文采用EPU算术平均值对前文计量方程重新进行回归，结果见表5列 （1）；采用综合创新指数替代创新产出水平对前文计量方程重新进行回归，结果见表5列 （2）。稳健性回归结果都显示经济政策不确定性与经济体创新产出显著负相关，与前文研究结果一致。
本文利用国内外现有三维生态足迹模型研究的理论成果，以生物生产型生态足迹理论为基础，通过改进三维生态足迹模型，对2010-2016年宁德市主要土地利用类型的生态足迹和生态承载力、生态足迹深度等进行分析，以期为宁德市经济、社会、生态可持续发展提供科学的理论参考和数据借鉴。
式中,S i 为景观类型i 的分离度,D i 为景观类型i 的距离指数,P i 为景观类型i 的面积指数,A i 为一个样区中景观类型i 的面积,n i 、A 含义同上。
(3)景观优势度指数DO i ：景观优势度表征景观结构中某一类型支配景观的程度,反映了该景观类型对景观格局形成和变化影响的大小[17]。公式为：
DO i =(1+景观类型i 的密度+景观类型i 的面积指数)/3
那么，我们应该如何选购对人体无害的樟脑丸呢？1.闻气味。天然樟脑丸有樟脑的特殊香味，而人工合成的樟脑丸则是刺鼻的味道。2.看成分表进行初步判断。天然樟脑丸会在成分中直接标明樟脑的含量，且价格较高；人工合成的樟脑丸则能在成分表中看到对二氯苯的含量，价格相对便宜。3.将樟脑球放入水中，可漂浮起来且不易融化的是天然樟脑丸，而人工合成的樟脑丸会沉于水底，且会逐渐融化。
将插值结果转化为栅格格式,并按照研究区范围进行不规则裁剪。通过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)。
2.3.2 景观脆弱度指数F i
近年来，“好教育”的推进使得广州教育的面貌发生了巨大转变。“十二五”时期，广州市顺利完成了“十二五”规划主要目标任务，广州好教育格局基本形成，市民群众对教育的获得感不断增强。《全国15个副省级城市教育现代化监测评价与比较研究报告（2015）》显示，4个教育现代化一级监测指标中，广州综合排名位居前列，其中教育普及发展第二、教育条件保障第二、教育质量要素第四。
景观脆弱度表示不同生态系统的易损性,与其在景观自然演替过程中所处的阶段有关。一般情况下,生态系统处于初级演替阶段时,食物链结构简单、生物多样性指数小,其较为脆弱。在土地利用生态系统中,人类活动是主要的干扰因素之一,土地利用程度不仅反映土地本身的自然属性,同时也反映人为因素与自然环境因素的综合效应。参考相关文献[11,17-18,23-26]并根据研究区实际情况,本研究区6种景观类型,以塌陷积水区最为脆弱,其次是其他农用地,最稳定的是城乡建设用地,6种景观类型分别赋值为：塌陷积水区为6、其他农用地为5、滩涂沼泽为4、耕地为3、水域较稳定为2、城乡建设用地为1。采用反正切函数归一化方法处理[11,17,23-25]之后,得到各景观类型的脆弱度指数F i 值分别为：0.8949、0.8743、0.8440、0.7952、0.7048、0.5。
2.3.3 土地利用生态风险指数ERI
基于上述所建立的景观干扰度指数和景观脆弱度指数,结合各景观类型的面积比重构建土地利用生态风险指数,用于描述一个评价单元内综合生态损失度的相对大小,以便通过采样方法将景观空间结构转化为空间化的生态风险变量。土地利用生态风险指数ERI计算公式如下：
式中,ERIk 为第k 个采样区内土地利用生态风险指数值,A ki 为第k 个采样区内土地利用类型i 面积,A k 为第k 个采样区面积,E i 为土地利用类型i 的干扰度指数,F i 为相应的脆弱度指数。
2.4 空间统计分析
2.4.1 空间自相关分析
空间统计学是以具有地理空间信息特性的事物或现象的空间相互作用及变化规律为研究对象,以区域化变量为基础,研究既具有随机性又具有结构性,或空间相关性和依赖性的一门科学[30-31],其核心是空间自相关性、空间依赖性和空间异质性。
所谓空间自相关性是指空间上越靠近的事物或现象越相似[4]。空间自相关的基本度量是空间自相关系数,可用全局和局部两种指标来度量。
全局空间自相关(Global Spatial Autocorrelation)用于描述区域单元某种属性值的整体分布状况,判断该属性值在空间上是否存在聚集性的特点。常用的指数是全局Moran′s I 指数。
全局Moran′sI 指数用来评估整个研究区域内所有空间对象是集聚分布、离散分布还是随机分布。其计算公式为[30-35]：
式中,n 为参与分析的空间单元数；x i 、x j 分别表示在空间单元i 和j 处空间对象的观测值
是空间权重矩阵。Moran′sI 的值域为[-1,1],大于0为正相关,小于0为负相关。数值越大表示空间分布的相关性越强,即空间上聚集分布的现象越明显；数值越小表示相关性弱,分散程度越高,即接近1表示聚集,接近-1表示离散；值趋于0时,代表空间分布呈随机分布现象。
2.4.2 地统计分析
(1)半变异函数
本研究借助地统计学中的半变异函数深入探讨景观结构变化等造成的土地利用生态风险的时空变化情况和空间变异规律。半变异函数的数学公式[4,32-36]为：
式中,r (h )是变异函数,h 是两样本配对抽样的分隔距离,即步长,Z (x i )和Z (x i +h )分别是随机变量Z 在空间位置x i 和x i +h 上的取值,N (h )为分隔距离为h 时的样本点总对数。
以r (h )为纵坐标,h 为横坐标做图,得到半方差图(semivariogram)。半方差图中包含3个重要参数：1)块金值(nugget)C 0,表示区域化变量在小于抽样尺度时的非连续性变异[35],2)基台值(sill)C +C 0,即平稳值,描述变量在研究区域范围中总的空间变异程度,基台值包括块金值(C 0)和结构方差(C )两部分。3)变程(range)A ,表示极限距离,即某特征在空间上自相关的空间幅度。
周教授一说，可蔓眼泪一下流了出来，说，真的穿越了哇？我好怕嘛。我要回去，我要回去找我妈嘛。说着就往外走。谷老板连忙拉住了，说，你瞎说个啥，你往哪儿走啊，你没看见外面啊？
(2)克里格插值
如果变异函数和相关分析的结果表明区域化变量存在空间相关性,就可以利用区域化变量的原始数据和变异函数的结构特点,对未采样点的区域化变量的取值进行线性无偏、最优估计[35]。
3 结果与分析
3.1 土地利用类型结构变化特征
分析不同尺度下多样性指数的半变异函数模型拟合参数(表4),尺度由0.5 km上升到1.5 km时,C 0从0.0228下降到0.0206,到2 km时C 0降至最低值0.0183。根据半变异函数的理论,块金值较大,说明较小尺度上的某种过程不可忽视,随机因素引起的空间异质性起主要作用。0.5 km到1 km C +C 0增加至0.0773,1 km到2 km C +C 0由0.0773降低至0.0633,说明随着尺度的增加,多样性指数的变化逐渐减弱。0.5 km到1 km C /(C +C 0)由0.606上升至0.730,说明0.5 km到1 km,多样性指数的空间自相关性逐渐增强；1 km到2 km C /(C +C 0),波动式下降,由结构性因素所引起的空间异质性在不断减弱。因此,1km能较直观的反映多样性指数的空间变异情况。
在防治玉米病虫害期间，可以利用化学除草的方式，在播种后可以进行封闭灭草，每亩用玉米草净100ml，并且兑水25kg，在播种之后的三天之内进行封闭灭草，另外土壤墒情较好时也可以采用播后苗前土壤封闭除草。
表1 2000—2016年研究区土地利用类型面积变化
Table 1 Area change of land use type in the research area from 2000 to 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′s I 指数分别为0.4576、0.4950、0.4279、0.4677,说明土地利用生态风险的空间分布呈集聚分布模式,具有很强的空间正相关性。
(1)景观破碎度指数C i ：景观破碎化是由于自然或人为干扰所导致的景观由单一、均质和连续的整体趋向于复杂、异质和不连续的斑块镶嵌体的过程,景观破碎度是景观异质性的重要组成部分,破碎度值越大,表明景观单元内部稳定性越低。公式为：
表2 1 km×1 km尺度下全局Moran′s I 指数及验证值
Table 2 Global Moran′s I index and test value at 1 km×1 km scale
3.3 土地利用生态风险的尺度效应分析
空间异质性依赖于尺度(粒度和幅度),粒度和幅度对空间异质性的测量和理解有重要的影响[4]。为了使建立的半变异函数模型能准确地反映各种尺度上的变化特征,要确定最佳的采样尺度[4,38]。根据土地利用生态风险指数的构成,以2005年土地利用类型分布图为试验数据,选取破碎度指数,并计算研究区Shannon多样性指数值,借助GS+地统计分析软件,进行不同尺度下的破碎度指数和Shannon多样性指数半变异函数模型拟合(表3、表4),以确定土地利用生态风险空间变异的最佳研究尺度。
表3 不同尺度下破碎度指数半变异函数拟合模型参数
Table 3 Fitting model parameters of semivariogram for fragmentation index at different scales
根据表3,0.5 km尺度下,破碎度指数的C 0为0.0014,C +C 0为0.0040,相比较1 km、1.5 km和2 km尺度下的C 0和C +C 0,其值最小,说明该幅度内尺度效应不明显；从0.5 km到1.5 km,C 0从0.0014上升到0.0056,说明人类活动等随机因素引起的空间异质性在逐渐增强[4],1.5 km到2 km C 0从0.0056骤然下降到0.0004,尺度的增加导致破碎度发生了明显变化。C +C 0从0.5 km到1 km上升最快,1 km到2 km缓慢上升,尺度的增加会降低空间自相关性。0.5 km到1 km C /(C +C 0)上升最快,之后缓慢上升,也反映了空间相关性在逐渐减弱。综上可知,0.5 km尺度过小,1.5 km到2 km空间自相关性明显减弱,1 km是适宜的研究尺度。
表4 不同尺度下Shannon多样性指数半变异函数拟合模型参数
Table 4 Fitting model parameters of semivariogram for Shannon diversity index at different scales
由表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；其他农用地和水域面积总体呈增加趋势,滩涂沼泽呈减少趋势。
综上,在进行微山县11个乡镇的土地利用生态风险时空异质性研究时,1 km是最佳研究尺度。
综上所述，知识可视化与网络媒介素养教育二者在互联网技术高度发达的信息社会中结合，能够实现网络技术更好地促进人的发展，实现人们在高技术环境下能够获得丰富的知识、方便的交流，从而实现自身的发展，最终实现智慧的技术环境，培养智慧的人的目标。
3.4 适宜尺度下土地利用生态风险的半方差分析
在GS+中,分别对4期1 km×1 km的样点数据进行半变异函数值的计算,并进行土地利用生态风险半变异函数模型拟合(表5及图4),其中步长设置为格网的间距1 km,有效滞后距离39 km为采样点最大间距的1/2。
其中,密度=(一个样区中景观类型i 的斑块数/此样区中总的斑块数),景观类型i 的面积指数同上。
表5 2000年—2016年土地利用生态风险半变异函数模型拟合参数
Table 5 Fitting model parameters of semivariogram for land use ecological risk from 2000 to 2016
由表5及图4看出,研究区2000年和2016年土地利用生态风险的最优拟合模型为球状模型,2005年和2010年的最优拟合模型为指数模型。4个时期的土地利用生态风险的空间变异特征主要表现在以下几方面：1)从各个年份最优拟合模型中块金值C 0的变化来看,2016年C 0最大,表明土地利用生态风险受人类活动等随机因素引起的空间异质性较大,较小尺度上的某种过程不能忽视。2)基台值C +C 0,从2000年的0.01000上升至2005年的最大值0.01276,到2016年有所降低,说明土地利用生态风险度空间变异程度随时间的推移在不断增强,其中2005年基台值相对突出,表明这个时段影响土地利用生态风险度的某些因子的空间变异性明显增强,2005年土地利用生态风险总的空间变异程度较高。3)结构方差与基台值的比值C /(C +C 0),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 插值预测误差统计表
Table 6 Interpolation prediction error statistics
3.6 土地利用生态风险的时空异质性分析
根据以上公式计算出C i 、S i 、DO i 指标后,由于量纲不同,需进行归一化处理。a 、b 、c 为各景观指数的权重,且三者相加为1。权重值的大小反映各景观指数解释景观类型受干扰的程度,根据研究区实际情况,并结合前人研究成果[17-18,23-24,26-29],认为破碎度指数最重要,其次为分离度指数和优势度指数,三个指数分别赋值为0.502、0.301、0.197。
图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 土地利用生态风险各年面积及比例
Table 7 Area and proportion of land use ecological risk from 2000 to 2016
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 土地利用生态风险各年面积变化及比例
Table 8 Change area and proportion of land use ecological risk from 2000 to 2016
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 土地利用生态风险时空面积变化及比例
Table 9 Spatial and temporal change area and proportion of land use ecological risk
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年较高,对此,在保证耕地总量动态平衡的前提下,应积极实行退耕还湖、退耕还湿的治理措施,增加芦苇沼泽湿地的面积,从而保护生物物种的多样性,提高生态系统的稳定程度。
中国水利：《意见》明确提出要建立水资源管理责任和考核制度，考虑到我国南北方水资源条件有较大差异，下一步如何进行责任分解和考核？
本研究基于景观生态学构建土地利用生态风险指数,但在景观干扰度指数的权重赋值时,没有进一步分析社会、经济因素与景观干扰度指数值间的定量关系,只是概念性的进行综合性度量,所以计算出的土地利用生态风险指数值也是相对的,并不具有绝对性；另外,在进行研究区土地利用生态风险半变异函数模型拟合时,虽然考虑了土地利用生态风险的尺度效应,但缺乏对生态风险的变化方向效应分析,未来可在这两个方面开展进一步的研究。
参考文献献(References):
[1] 卢宏玮, 曾光明, 谢更新, 张硕辅, 黄国和, 金相灿, 刘鸿亮. 洞庭湖流域区域生态风险评价. 生态学报, 2003, 23(12): 2520-2530.
[2] U.S. Environmental Protection Agency. Framework for Ecological Risk Assessment. Washington, DC: Risk Assessment Forum, 1992: 1-2.
[3] 付在毅, 许学工. 区域生态风险评价. 地球科学进展, 2001, 16(2): 267-271.
[4] 邬建国. 景观生态学: 格局、过程、尺度与等级(第二版). 北京: 高等教育出版社, 2007: 19-19, 126-136.
[5] 付在毅, 许学工, 林辉平, 王宪礼. 辽河三角洲湿地区域生态风险评价. 生态学报, 2001, 21(3): 365-373.
[6] 许学工, 林辉平, 付在毅, 布仁仓. 黄河三角洲湿地区域生态风险评价. 北京大学学报: 自然科学版, 2001, 37(1): 111-120.
[7] 王春梅, 王金达, 刘景双, 张学林, 严登华. 东北地区森林资源生态风险评价研究. 应用生态学报, 2003, 14(6): 863-866.
[8] 程建龙, 陆兆华, 范英宏. 露天煤矿区生态风险评价方法. 生态学报, 2004, 24(12): 2945-2950.
[9] 陈辉, 李双成, 郑度. 基于人工神经网络的青藏公路铁路沿线生态系统风险研究. 北京大学学报: 自然科学版, 2005, 41(4): 586-593.
[10] Breyta R, Black R, Kaufman J, Kurath G. Spatial and temporal heterogeneity of infectious hematopoietic necrosis virus in Pacific Northwest salmonids. Infection, Genetics and Evolution, 2016, 45: 347-358.
[11] 李谢辉, 王磊, 李景宜. 基于GIS的渭河下游河流沿线区域生态风险评价. 生态学报, 2009, 29(10): 5523-5534.
[12] 常青, 邱瑶, 谢苗苗, 彭建. 基于土地破坏的矿区生态风险评价: 理论与方法. 生态学报, 2012, 32(16): 5164-5174.
[13] 许妍, 高俊峰, 郭建科. 太湖流域生态风险评价. 生态学报, 2013, 33(9): 2896-2906.
[14] 张雅洲, 谢小平. 基于RS和GIS的南四湖生态风险评价. 生态学报, 2015, 35(5): 1371-1377.
[15] 肖武, 李素萃, 王铮, 杨耀淇, 王涛. 高潜水位煤矿区生态风险识别与评价. 生态学报, 2016, 36(17): 5611-5619.
[16] 巫丽芸, 何东进, 游巍斌, 邓西鹏, 谭勇, 纪志荣. 福建东山岛灾害生态风险的时空演化. 生态学报, 2016, 36(16): 5027-5037.
[17] 吴健生, 乔娜, 彭建, 黄秀兰, 刘建政, 潘雅婧. 露天矿区景观生态风险空间分异. 生态学报, 2013, 33(12): 3816-3824.
[18] 李保杰. 矿区土地景观格局演变及其生态效应研究——以徐州市贾汪矿区为例[D]. 徐州: 中国矿业大学, 2014.
[19] 胡金龙, 周志翔, 滕明君, 罗楠. 基于土地利用变化的典型喀斯特流域生态风险评估——以漓江流域为例. 应用生态学报, 2017, 28(6): 2003-2012.
[20] 孙洪波, 杨桂山, 苏伟忠, 朱天明, 万荣荣. 沿江地区土地利用生态风险评价——以长江三角洲南京地区为例. 生态学报, 2010, 30(20): 5616-5625.
[21] 王国杰, 廖善刚. 土地利用强度变化的空间异质性研究. 应用生态学报, 2006, 17(4): 611-614.
[22] 曾勇. 区域生态风险评价——以呼和浩特市区为例. 生态学报, 2010, 30(3): 668-673.
[23] 谢花林. 基于景观结构和空间统计学的区域生态风险分析. 生态学报, 2008, 28(10): 5020-5026.
[24] 谢花林. 基于景观结构的土地利用生态风险空间特征分析——以江西兴国县为例. 中国环境科学, 2011, 31(4): 688-695.
[25] 肖琳, 田光进. 天津市土地利用生态风险评价. 生态学杂志, 2014, 33(2): 469-476.
[26] 赵丽红. 南昌市景观格局时空变化及其驱动力研究[D]. 南昌: 江西农业大学, 2016.
[27] 于化龙, 陈青锋, 田超, 张杰, 陈亚恒. 基于景观结构的土地利用生态空间特征及风险评估——以怀来县为例. 水土保持研究, 2016, 23(3): 155-163.
[28] 卿凤婷, 彭羽. 基于景观结构的北京市顺义区生态风险时空特征. 应用生态学报, 2016, 27(5): 1585-1593.
[29] 王鹏, 王亚娟, 刘小鹏, 陈晓, 孔福星. 基于景观结构的生态移民安置区生态风险评价——以宁夏红寺堡区为例. 生态学报, 2018, 38(8): 2672-2682.
[30] 胡明星, 李建. 空间信息技术在城镇体系规划中的应用研究. 南京: 东南大学出版社, 2009: 80-80, 82-94.
[31] 冯益明. 空间统计学理论及其在林业中的应用. 北京: 中国林业出版社, 2008: 1-4.
[32] 闫庆武. 空间数据分析方法在人口数据空间化中的应用. 南京: 东南大学出版社, 2011: 68-74, 97-114.
[33] 刘湘南, 黄方, 王平. GIS空间分析原理与方法(第二版). 北京: 科学出版社, 2008: 189-206.
[34] 吴剑, 陈鹏, 文超祥, 傅世锋, 陈庆辉. 基于探索性空间数据分析的海坛岛土地利用生态风险评价. 应用生态学报, 2014, 25(7): 2056-2062.
[35] 徐建华. 计量地理学(第二版). 北京: 高等教育出版社, 2014: 158-181.
[36] 李哈滨, 王政权, 王庆成. 空间异质性定量研究理论与方法. 应用生态学报, 1998, 9(6): 651-657.
[37] 邓敏, 樊子德, 刘启亮. 空间分析实验教程. 北京: 测绘出版社, 2015: 66-70.
[38] 岳文泽, 徐建华, 徐丽华, 谈文琦, 梅安新. 不同尺度下城市景观综合指数的空间变异特征研究. 应用生态学报, 2005, 16(11): 2053-2059.
Spatial-temporal heterogeneity of ecological risk of land use in mining areas
CHANG Xiaoyan1,2, LI Xinju 2, *, LI Xican1, GUO Peng1, GAO Feng3
1College of Information Science and Engineering ,Shandong Agricultural University ,Tai ′an 271018,China 2College of Resources and Environment ,Shandong Agricultural University ,Tai ′an 271018,China 3Ji ′ning Land Resources Bureau ,Shandong Province ,Ji ′ning 272001,China
Abstract ： In this study, the 11 towns with the most coal mines and serious surface collapses in Weishan County were selected as the areas of study, and remote sensor images and land use data were used as the data sources. From this information, distribution maps of land use type in 2000, 2005, 2010 and 2016 were obtained. Based on landscape ecology theory, an index of the ecological risk of land use was constructed. In this study, the theories of spatial statistics and geostatistics were used to determine an appropriate study scale. Spatial-temporal change and spatial heterogeneity of the ecological risk of land use in the research area were analyzed between 2000 and 2016. It was concluded that 1 km × 1 km regular grid was the best research scale. At this scale, the spatial distribution of the ecological risk of land use in the research area was concentrated and had a strong positive spatial correlation. The random variation of the land use ecological risk index on a small scale was gradually replaced by spatial structural variation on a larger scale between 2000 and 2005. The total spatial variability of ecological risk in 2005 was high. Between 2005 and 2016, the degree of spatial variability caused by random factors was increased. Between 2000 and 2016, the ecological risk of land use was characterized by the transition from medium-high and high ecological risk to medium and medium-low ecological risk. The ecological risk of deteriorated areas and obviously deteriorated areas was mainly distributed around the eastern mining area, roadside, and serious collapse areas. In addition, the ecological risk of deterioration of cultivated land on both sides of the east and west side of Weishan Lake was higher between 2010 and 2016. Over the past 17 years, the stable areas and improved areas of Weishan Lake increased significantly, with a clear decrease in the overall ecological risk.
Key Words : land use; landscape pattern; ecological risk; scale; spatial autocorrelation; spatial-temporal heterogeneity
DOI: 10.5846/stxb201803290636
常小燕,李新举,李西灿,郭鹏,高峰.矿区土地利用生态风险的时空异质性.生态学报,2019,39(9):3075-3088.
Chang X Y, Li X J, Li X C, Guo P, Gao F.Spatial-temporal heterogeneity of ecological risk of land use in mining areas.Acta Ecologica Sinica,2019,39(9):3075-3088.
基金项目: 山东省重点研发计划(2016ZDJS11A02); 国家自然科学基金项目(41771324)
收稿日期: 2018-03-29; 网络出版日期:2019-02-27
*通讯作者 Corresponding author.E-mail: lxj0911@126.com
标签：土地利用论文; 景观格局论文; 生态风险论文; 尺度论文; 空间自相关论文; 时空异质性论文; 山东农业大学信息科学与工程学院论文; 山东农业大学资源与环境学院论文; 山东省济宁市国土资源局论文;