生态安全是21世纪人类社会可持续发展所面临的一个新主题^[1], 土地生态安全作为生态安全问题研究的重要组成部分, 己成为当前土地资源可持续利用研究的前沿课题。近年来, 由于我国人口急剧增加、工业化和城镇化进程加快, 人地矛盾愈发凸显, 不合理的土地利用方式以及土地利用强度的不断加大, 导致部分区域土地生态系统遭受到严重的威胁和破坏, 土地生态安全问题成为区域土地资源可持续利用的主要“瓶颈”。因此, 如何保障区域土地生态系统的安全、实现土地生态安全预警成为关系到人类社会能否持续稳定发展的关键问题。

土地生态安全预警研究起源于生态安全预警, 可以追溯到20世纪70年代^[2]。随着1975全球环境监测系统(GEMS)的建立^[2], 生态安全预警研究开始得到广泛关注；在此基础上部分学者^[3-4]从不同角度对生态环境预警进行了开拓性研究, 生态环境预警理论不断完善, 技术方法和手段不断更新和提高, 从单项预警发展到综合预警, 从专题预警发展到区域预警^[5]。随着生态环境预警理论研究的不断深入, 生态环境预警理论被逐步引入到土地科学领域, 土地生态安全预警随着土地生态安全问题的产生和信息技术的发展亦随之出现^[6], 国内外学者进行了大量研究。在国外土地生态安全预警方面, 相关研究以监测预警为主, 主要包括土壤监测系统建设^[7], 土地利用变化广域性、综合性的动态监测^[8-11]和农业资源环境、土地生态环境预警预测等方面^[12-15]。在国内土地生态安全预警方面, 主要涉及土地生态安全预警理论^{[6, 16-18]}、预警指标体系^[19-22]、预警方法^[23-26]及预警信息系统建设^[27-28]等方面。综合来看, 尽管国内外学者针对土地生态安全预警开展了大量研究, 但仍存在以下不足：一是土地生态安全预警多以县(区)行政单位为评价单元, 难以反映县(区)内部土地生态安全状况的局部差异；二是土地生态预警评价指标体系过多地依赖统计年鉴等社会统计资料, 对土地生态状况的自然属性重视不够；三是忽略了土地生态安全警度在空间上的聚集规律、关联模式等空间分布特征及其演变过程, 而土地生态安全警度的时空演变规律更具实用性, 更有利于揭示区域土地生态状况的变化特征及驱动机制, 能更好地为土地生态环境保护及人地关系协调发展提供参考。

三峡库区生态敏感区位于长江流域生态屏障的咽喉地带, 是中国具有全球保护意义的生物多样性关键地区之一, 同时也是长江上游经济带西部中心枢纽, 保护资源与保障发展的矛盾日益突出。本研究选取三峡库区生态敏感区典型区域万州区2000年、2009年、2014年3个时期以行政村为评价单元开展土地生态安全预警测度与时空演变研究, 构建基于PSR模型的区域土地生态安全预警指标体系, 应用空间自相关模型和GIS相关技术开展区域土地生态安全警度时空分异特征, 旨在为破解三峡库区生态敏感区生态环境瓶颈要素制约, 实现区域土地资源可持续利用和社会经济可持续发展提供决策参考。

1 研究区概况

重庆市万州区位于三峡库区腹心地带、四川盆地东部边缘, 地处107°52′22″—108°53′25″E, 30°24′25″—31°14′58″, 东与云阳, 南与石柱和湖北省利川, 西与忠县和梁平, 北与开县和四川省达州接壤。全区辖11个街道办事处, 29个镇, 12个乡, 幅员面积3456.55 km²。万州区在地势上处于我国第二阶段的东部边缘地带, 四川盆地向秦巴山地和云贵高原的过渡带, 境内地势总体上以长江为界, 长江以北为北高南低, 长江以南为南高北低, 区内山脉平行延伸, 谷岭相间分布, 山脊线构造线基本一致, 多呈东北—西南走向, 形成以中低山和丘陵为主的地貌类型, 属亚热带季风湿润带, 气候温和, 四季分明, 冬暖夏热, 无霜期长, 雨量充沛, 年均降雨量1416 mm。万州区素有“川东门户”之称, 是辐射川东、渝东、鄂西、黔东北、陕南、湘西的水陆交通枢纽和物资集散中心, 区位优势十分明显。随着国家长江经济带战略的深入实施, 万州区作为三峡库区腹心地带的区域性中心城市, 其土地生态安全建设不仅是区域经济可持续发展的目标, 更是三峡库区国土资源安全和生态文明建设的重要保障。

2 数据来源与处理方法 2.1 数据来源

本研究土地利用现状数据以万州区2000年、2009年、2014年土地利用矢量图、1:1万DEM、土地利用现状图、行政区划图为基础, 同时结合区域土壤、气候、植被等自然地理资料, 利用ArcGIS 10.2、GeoDa等软件对土地利用数据进行分析处理。其他社会经济数据来源于2001年、2010年、2015年万州区统计年鉴及林业、环保、农业等部门统计资料和实地调查资料。

2.2 土地生态安全预警评价模型 2.2.1 预警指标体系构建

本研究遵循科学性、综合动态性、指标可获得及可操作性、灵敏性原则, 综合三峡库区生态敏感区土地生态状况实际情况^[29-31], 参照20世纪80年代末联合国环境规划署(UNEP)和经济合作开发组织(OECD)联合提出的“压力一状态一响应”(Pressure-State-Response, 即PSR)框架模型和联合国(UN)修改PSR模型后提出的“驱动力一状态一响应”(Driving Force-State-Response, 即DSR)框架模型^[32-34], 并结合《生态环境状况评价技术规范》(HJ 192—2015), 采用德尔菲法、相关分析法, 充分运用国土、环保、林业、农业等领域15位专家深厚的专业理论知识和丰富的实践经验, 选取了21项指标构建了预警指标体系(表 1)。本研究构建的PSR土地生态安全预警指标体系中, “压力指标”是指人类活动对土地生态系统安全产生的压力, 具体由人口、经济、城镇化及交通建设发展等方面的压力构成；“状态指标”用来反映特定时间阶段土地生态安全状况与动态特征, 包括土地生态系统的自然环境状态以及土地利用状况等方面；“响应指标”是能够反映处理土地生态安全问题和维护改善土地生态安全状态的能力等方面。

2.2.2 指标阈值确定与标准化

参照相关研究^{[26, 35]}, 将三峡库区生态敏感区土地生态安全等级划分为无警、轻警、中警、重警、巨警五级。在确定土地生态安全预警指标体系基础上, 参考国家、行业和地方规定的相关环境质量标准, 结合相关研究成果^{[18, 30, 36-37]}, 从三峡库区土地生态环境现状与特征^[29-31]出发, 采用绝对确定法、相对确定法、专家确定法、借鉴法^[37]和自然分类法(该方法利用统计学Jenk最优化法得出分界点, 能够使各级内部方差之和最小)^[38]确定各预警指标对应的各生态安全级别的阈值(表 1), 并据此阈值采用内插法进行预警指标标准化处理。

2.2.3 权重确定

本研究选取熵权法计算预警指标权重, 并选取国土、环保、林业、农业等领域15位专家采用特尔菲法进行校核, 综合确定区域土地生态安全预警指标权重(表 1), 熵权法权重计算过程如下：

1) 计算第j项指标的信息熵e_j:

(1)

其中, k:调节系数, , 且k>0；P_ij:第i个参评对象第j个评价指标的标准化值。

(2) 计算第j项指标的效应值g_j:

(2)

(3) 计算第j项指标的权重w_j:

(3)

某项评价指标的熵值越大, 其权重值越小, 反之亦然, 研究区土地生态安全预警指标的权重值详见表 1。

2.2.4 土地生态安全预警指数测算

参照相关研究^[39-40], 采用乘算模型对评价单元进行土地生态安全预警指数计算, 公式如下：

(4)

式中：ESWI：某类指标的土地生态安全预警指数；A_i为评价指标标准化值；W_i为土地生态安全预警指标A_i的权重；n为指标总项数。

2.2.5 警度判别标准

本研究采用系统聚类法计算模糊等价关系进行初始分类, 再经ISODATA聚类法反复迭代进行研究区土地生态安全警度分类, 构造无警点、低警点、中警点、重警点、巨警点等多个警限, 并划分出无警、轻警、中警、重警、巨警等5个土地生态安全预警等级, 详见表 2。

2.3 土地生态安全警度时空演变统计学方法 2.3.1 地统计学方法

目前, 地统计学方法已成为研究生态安全空间变异性的一种有效方法。区域土地生态安全预警指数是一种空间变量, 空间变化特征具有结构性和随机性, 可以采用地统计学中空间变异函数的方法, 借助半方差函数进行区域土地生态安全警度的空间分析^[41]计算公式如下：

(5)

式中：r(h)是半变异函数, 揭示了整个尺度上的空间变异格局；h为两样本点的空间分割距离；Zx_i和Zx_i+h分别为空间位置x_i和x_i+h处的数值；n(h)为空间距离为h时的样本对数。

2.3.2 空间自相关

地理学第一定律指出：任何事物之间均相关, 而离得较近的事物总比离得较远的事物相关性要高^[42], 亦即地理实体由于受到空间吸引和空间扩散的影响而表现出一定的相关性^[43]。空间自相关分析是理论地理学的基本方法之一^[44], 主要用来分析空间数据的统计性分布规律。根据空间分析尺度不同, 土地生态安全空间自相关可以划分为全局型空间自相关和局域型空间自相关两种。全局指标用于验证整个研究区域某一要素的空间相关关系, 而局部指标则用于反映整个大区域中的一个局部小区域单元上的某种地理现象或某一属性与相邻局部小区域单元上同一现象或属性的相关程度。本研究运用全局空间自相关Moran′s I指标和局部空间自相关指标LISA(Local Indicators of Spatial Association)来分析万州区土地生态安全警度的空间分布特征。Moran′s I和LISA指标的计算公式如下^[45-47]：

(6)

式中, , , x_i表示第i地区的观测值, n为栅格数, W_ij为二进制的邻接空间权重矩阵, 表示空间对象的邻接关系。i=1, 2, …, n；j=1, 2, …, m；当区域i和区域j相邻时, W_ij=1；当区域i和区域j不邻时, W_ij=0。Moran′s I一般介于-1到1之间, 大于0为正相关, 小于0为负相关, 绝对值越大表示空间分布具有越强关联性, 即空间上有强聚集性或强相异性。反之, 绝对值越小, 表示空间分布关联性较小, 当值趋于0时, 即代表此时空间分布呈随机性。

局部空间自相关Local Moran′s I(LISA)是将Moran′s I分解到各个空间单元, 其公式为：

(7)

式中, , I_i值为正时表示该空间单元周围相似值(高值或低值)的空间集聚, I_i值为负时表示非相似值之间的空间集聚。

3 结果统计与分析 3.1 区域土地生态安全警度时间演变特征

由图 1和表 3可以看出, 2000年土地生态安全警度以无警和轻警为主, 两者之和占到了区域土地总面积的78.1%, 说明了区域土地生态状况总体情况较好；而警度较高区域(重警和巨警)主要集中分布在中心城区所在乡镇及其邻近区域(周家坝街道、百安坝街道、牌楼街道、沙河街道、太白街道、高笋塘街道、钟鼓楼街道、陈家坝街道、龙都街道、五桥街道、九池乡、双河口街道), 仅占区域土地总面积的7.08%。2009年区域土地生态安全发生了较大变化, 其土地生态安全警度明显提高, 警度较低区域(无警、轻警)面积比例下降到48.34%, 而警度较高区域(重警和巨警)面积比例则上升到18.47%, 除了2000年土地生态安全警度较高区域外, 大周镇、小周镇、天城镇、太龙镇、高梁镇、熊家镇等中心城区邻近区域亦由2000年的低警度区域变为2009年的高警度区域。2014年与2009年相比, 土地生态安全警度变化不明显, 略有提升, 警度较低区域(无警、轻警)面积比例进一步下降到45.61%；而警度较高区域(重警和巨警)面积比例则由2005年的18.47%下降到2014年的15.63%, 这是由于部分高警度区域(高梁镇、熊家镇、太龙镇等)土地生态环境质量得到一定程度提升所致, 其生态安全警度由2009年的重警转变为2014年的中警状态。总体来说, 研究区2000—2014年土地生态安全警度呈连续上升趋势, 警度较高区域由万州区中心城区不断向四周扩散, 但2000—2009年期间研究区土地生态安全警度提升幅度比2009—2014年期间更为明显, 可见2000—2009年期间是区域土地生态安全警度变化的关键时段。这主要是由于2000—2009年期间区域城镇化和工业化进程的加快, 加剧了人地矛盾, 农田、园地和未利用地都存在被不同程度开发占用的现象^[48], 进而导致土地生态环境受到越来越大的压力, 生态系统稳定性逐步变差, 土地生态安全警度亦随之逐步提升。而随着近年来国家生态文明建设战略的大力实施和绿色发展理念的深入践行, 研究区土地生态状况得到的保护力度不断增强, 这在一定程度上缓解了区域2009—2014年期间土地生态安全警度快速提升的趋势。

3.2 区域土地生态安全警度空间结构分析 3.2.1 空间结构变异性分析

结合土地生态安全预警指数, 运用GS+地统计学软件对研究区土地生态安全警度进行变异函数理论模型的拟合, 模型拟合结果见表 4。可以看出, 2000年、2009年、2014年研究区土地生态安全预警指数均为球形模型的拟合效果最佳, 复相关系数分别R²分别为0.838、0.949、0.961。

土地生态状况空间分布具有异质性, 其空间分布受结构性因素和随机性因素(或称非结构性因素)的双重作用, 可通过拟合空间结构模型来探索区域土地生态安全警度空间格局变化的内在机理。块金值C₀表示距离为零时的空间变异程度, 反映了区域化土地生态状况的随机性特征, 块金值越大则区域土地生态状况分布随机性越强。基台值(C₀+C)表示空间的最大变异程度, 基台值越高, 区域土地生态状况的空间异质性越强。C₀/(C₀+C)表示随机性变异占总变异的比例, 该比值越高, 随机性因素影响越大。随机性因素(外因)包括各种自然灾害以及人为活动等因素；而结构性因素(内因)包括气候、地形地貌、土壤类型、植被类型等主导区域土地生态安全状况的空间分布因素。变程A₀反映了区域土地生态状况空间自相关尺度的范围, 若要素间距离在A₀范围内, 则要素间存在空间自相关；否则, 要素呈随机性分布。

2000—2014年研究区土地生态安全警度的块金效应逐步增强, 从2000年的0.00061增加到2014年的0.001, 表明研究区随机性因素逐步增强, 不可忽视。C₀/(C₀+C)在2000年、2009年、2014年分别为13.8%、18.07%、19.34%, 表明随机性因素对区域土地生态安全警度空间分异影响效应逐步增强, 这与研究区14年来经济高速发展、人类对土地资源利用强度加剧以及区域自然灾害频发的实际情况相符。2000年、2009年、2014年研究区土地生态安全警度空间分异变程分别为2.26、3.43、3.59 km, 可以看出研究区土地生态安全警度的相关性范围亦在不断扩张。

3.2.2 空间结构关联性分析

(1) 全局自相关分析

利用Moran′s I指数判定区域土地生态安全预警指数的空间相关性和聚集程度。2000年、2009年、2014年研究区域土地生态安全预警指数的Moran′s I值分别为0.7823、0.7772、0.7750, 整体呈略微下降趋势, 表明研究区域土地生态安全警度在空间上具有较强的正相关。通过Moran散点图可以定性区分某一区域与周边区域某种属性值的相互关系。Moran散点图分为4个象限, 第一象限表示高—高聚集(HH)、第二象限表示低—高聚集(LH)、第三象限表示低—低聚集(LL)、第四象限表示高—低聚集(HL)。由图 2可知, 三个年份土地生态安全警度主要集中在HH和LL区, 少量分布于HL、LH两个区, 表明研究区土地生态安全警度具有较强的空间正相关性, 空间集聚特征显著, 亦即土地生态安全警度在空间上存在趋于集群的现象。

(2) 局部自相关分析

全局Moran′s I指数可检验研究区土地生态安全警度的总体空间分布模式, 却不能反映相邻区域间的要素或属性的空间关联模式及局域显著性水平, 而局域空间自相关能更深入地发现集聚和变异特征^[49], 因此需要通过局域指标来进一步研究相邻区域之间的土地生态安全警度相关程度, 进而深入探讨区域土地生态安全警度的空间格局及其可能成因。通过对研究单元的土地生态安全预警指数进行局域空间关联分析, 得到LISA集聚图(图 3)、LISA显著性检验图(图 4)和LISA面积统计表(表 5)。

从图 3和表 5可以看出, 2000年“高—高”值分布集中, 主要集中分布在中心城区所在乡镇及其邻近区域(周家坝街道、百安坝街道、牌楼街道、沙河街道、太白街道、高笋塘街道、钟鼓楼街道、陈家坝街道、龙都街道、五桥街道、九池乡、双河口街道、高峰镇、高粱镇、天城镇、熊家镇、长岭镇等), 大周镇、小周镇、太龙镇、新乡镇、走马镇、地保乡略有零星分布。中心城镇及其邻近区域历来是万州区经济社会快速发展地区, 地势平坦、土壤质量好, 人类土地利用活动剧烈, 土地生态安全警度较高；“低—低”值分布相对零散, 主要分布在研究区西部的武陵镇、分水镇、孙家镇、柱山乡、弹子镇和东南部的新田镇、茨竹乡、龙驹镇、恒合乡, 与这一时期的土地生态安全警度较低区域(无警、轻警)分布基本一致；“低—高”值主要集中分布在钟鼓楼街道、天城镇、熊家镇、高粱镇、小周镇等中心城区所在乡镇及其邻近区域；“高—低”值仅有零星分布。

2000—2009年研究区土地生态安全警度“高—高”值地区基本没有发生变化, 但是范围进一步扩展, 主要是大周镇、小周镇、太龙镇等乡镇的部分不显著区域和熊家镇、天城镇、高粱镇等乡镇的“低—高”值区变为了2009年的“高—高”值区。“低—低”值总体分布情况与2000年基本一致, 主要分布在万州区西部的弹子镇、余家镇、孙家镇、后山镇、分水镇、响水镇、郭村镇、柱山乡和万州区东南部的茨竹乡、龙驹镇、恒合乡、罗田镇等地区。

2009—2014年研究区土地生态安全警度“高—高”值范围较2009年有所缩小, 主要体现是高粱镇、太龙镇等部分“高—高”值区变为了不显著区域和熊家镇等部分“高—高”值区变为了“低—高”值区；“低—低”值区未发生明显变化, 仍主要集中分布在研究区西北部和东南部等土地利用程度相对较弱、植被覆盖相对较好地区。可以看出, 研究区“高—高”值区总体上呈现出2000—2009年范围逐渐增加, 而2009—2014年范围缩小的趋势, 这与区域土地生态安全警度变化情况相一致。土地生态安全警度的空间分布并不是随机的, 存在空间上的聚集趋势；土地生态安全警度较低的单元可以对周边地区产生正面影响, 使其周围土地生态警度降低, 而“高—高”值区也可以对周边地区产生负面影响, 使其土地生态安全警度提高。

从图 4和表 5可以看出, 土地生态安全警度集群分布的显著性水平分异明显, “高—高”值、“低—低”值区均以极显著水平(P < 0.01)为主, 而“低—高”值区、“高—低”值区则均以显著水平(P < 0.05)为主。研究区域2000—2014年土地生态安全警度集群分布的显著性水平发生了明显改变, 其“高—高”值区极显著性水平的区域面积从2000年的36934.24 hm²增加到2009年的43132.95 hm², 再减少到2014年的40596.52 hm², 呈现出先增后减、总体增加的趋势, 与区域土地生态安全警度2000—2014年期间的总体提升, 而2000—2014年期间略有缓解的总体变化特征相一致。

4 结论与讨论

本研究以三峡库区生态敏感区典型区域万州区为研究对象, 基于PSR模型构建土地生态安全预警指标体系, 计算土地生态安全预警指数并划分出土地生态安全预警等级, 通过半变异函数模型和空间自相关分析方法, 以2000年、2009年、2014年为3个年份为基础, 分析2000—2014年研究区土地生态状况警度时空演变状况, 结论及建议如下：

(1) 本研究选取的土地生态安全预警指标是基于三峡库区生态敏感区自然生态基底与人类社会发展对生态环境需求而选取的评价指标, 是能较好地反映研究区土地生态状况实际的个性化指标, 改变了以往区域土地生态安全评价中过多地依赖于社会经济统计资料而对土地生态状况自然属性重视不够的不足之处, 以更加科学合理地表征研究区土地生态安全警度。

(2) 本研究以土地利用现状数据库中行政村为基本评价单元, 能很好地展示三峡库区生态敏感区典型区域万州区的土地生态安全警度的时空分异状况, 极大地提高了区域土地生态状况评估精度, 可以为整个三峡库区生态敏感区土地生态安全时空分异规律提供参考依据。研究表明, 研究区2000—2014年土地生态安全状况呈现连续下降趋势, 土地生态安全警度较高区域由万州区中心城区不断向四周扩散, 但2000—2009年期间研究区土地生态安全状况变化比2009—2014年期间明显, 这主要得益于近年来国家生态文明建设战略的大力实施和绿色发展理念的深入践行, 区域土地生态状况得到的保护力度不断增强, 这在一定程度上缓解了2009—2014年期间土地生态状况快速下降的趋势。可以看出, 尽管区域土地生态状况快速恶化的趋势得到一定程度缓解, 但仍不容乐观, 为改善区域土地生态状况, 在“十三五”期间乃至今后更远时期内, 有关部门必须在保护耕地、控制建设用地盲目扩张和水土流失、预防地质灾害、减少污染物排放、增加绿化建设、提高绿色发展水平等方面做出更加积极有效的排警调控措施, 促进土地生态系统安全平稳运行, 为社会经济可持续发展提供坚实保障。

(3) 研究表明, 土地生态安全警度在空间上存在趋于集群的现象, “高—高”值、“低—低”值区均呈现出2000—2009年范围逐渐增加, 而2009—2014年范围逐渐缩小的趋势, 且“高—高”值区主要集中分布在中心城区所在乡镇及其邻近区域, “低—低”值区则主要集中分布在研究区西北部和东南部等土地利用程度相对较弱、植被覆盖好地区。近年来, 随着重庆市“五大功能区”发展战略的深入实施和新型工业化、城镇化战略的快速推进, 区域土地利用结构发生了显著变化, 城市周边土地生态安全警度的聚集趋势会愈发明显。因此, 中心城区周边的“低—高”值区将是城市发展中优先被占用的区域, 在未来的土地利用中, 尤其需加强对城区周边“低—高”值区的重点保护。

(4) 构建基于PSR模型的土地生态安全预警指标体系和采用地统计学方法、空间自相关分析方法的研究适用于区域土地生态安全警度变化的诊断, 有利于促进区域土地利用生态环境的保护和改善。本研究对进一步完善区域土地生态安全预警具有一定的参考价值, 对于当前受到高度关注的国家生态安全战略问题的认知与解决具有重要意义, 同时也为区域土地资源利用的规划、管理以及社会经济的可持续发展提供决策参考。但是, 区域土地生态安全预警是一项系统工程, 本研究设置的土地生态安全预警指标阈值等问题是一个长期复杂的学术问题, 目前学术界对这个问题还没形成一定的定论, 有待深入研究。