土地是人类赖以生存的基础，土地生态安全与人类的生存生产紧密相关^[1]。土地生态安全是指在一定时空范围内，土地生态系统能维持其结构性与功能性健康的前提下，满足人类生存和发展所必需的资源供给和生态需要的能力^[2-3]。土地生态安全预警作为土地生态安全研究的重要分支，是对土地生态环境变化以及土地生态环境与社会经济协调变化的警告，预前揭示潜在的风险或危害，对维护区域土地生态安全具有重要意义^[4]。一直以来的高强度土地开发与利用，极大的破坏了土地生态平衡，土地生态安全形势愈发严峻。现有研究从系统分解综合考量生态系统与社会经济系统间相互关系为起点，构建适宜土地生态安全预警指标体系，对研究单元逐项评价后经空间叠加、模型计算等实现对区域土地生态安全的相对测度^[5-6]。建立科学合理的预警指标体系是土地生态安全定量测度的基础，也是土地生态安全预警的关键^[7]。目前，构建土地生态安全预警框架体系主要有3种：PSR(压力—状态—响应)概念模型、EES(经济—环境—社会)概念模型以及基于生态系统本底特征指标体系^[8]。PSR概念模型的衍生模型(如PISR，DPSIR，DSR和DPSER等)、暴露—响应概念性框架模型、基于土地子系统—耕地的生态安全指标体系等^[6]。研究尺度主要包括行政区、城市(圈/群)^[9-10]、重点流域^[11-12]、特殊的生态敏感区等^[13-14]。综合来看，现有成果为后期研究提供了科学范式，但土地生态安全存在明显的空间相关性及异质性特征，仅仅利用传统计量模型进行影响因素分析，无法揭示其地理空间时空规律，而运用空间自相关模型及地理探测器可以弥补该问题的不足，充分揭示土地生态安全问题的驱动机制及发生原理。

汾河是黄河流域第二大流域单元，也是山西省重要的生产生活生态空间集聚区，其占有全省约25%土地面积、26%的水资源，却承载约40.1%的人口和42.6%的GDP^[15]。长期以来高强度的人类扰动造成流域形成高达2.48×10⁴ km²水土流失区，同时还存在土壤污染、植被破坏等一系列的土地生态安全问题^[16]。然而，目前形成该问题的驱动机制、发生原理尚不明确，所以，本文从土地生态安全时空格局演变角度，以DPSIR理论框架构建研究区预警指标体系，采用乘算模型定量测算2006，2013，2018年基准期的安全预警指数的基础上，运用空间自相关模型及地理探测器揭示空间相关性及异质性机理，探明研究区土地生态安全警度时空演变的驱动机制及发生原理，为流域土地资源科学合理利用及可持续发展提供参考与依据。

1 研究区概况及研究方法 1.1 研究区概况

汾河源于宁武县管涔山，在万荣县注入黄河，全长约713 km。汾河流域地处山西省中部(35°22′6″—38°49′24″N，110°30′3″—112°50′24″E)，黄土高原东缘，地势由北向南倾斜，海拔高度240~2 786 m，地貌以盆地、平川、丘陵为主，并由地堑型纵谷串联太原、临汾等盆地构成。属温带大陆性季风气候，年均气温11 ℃，年均降水量392.8 mm。土壤类型含潮土、褐土、棕壤、粗骨土、红黏土等8类，褐土为主要土壤之一。流域总面积约3.95万km²，涉及忻州、太原、晋中、吕梁、运城、临汾6市40县/区，目前，流域集中全省约40.1%的人口，42.6%的GDP，城镇化水平达60.3%，为全省重要的生产、生活和生态空间聚集地。

1.2 数据来源

本文社会经济统计等数据来源于2006, 2013, 2018年《山西省统计年鉴》^[17]《山西省水资源公报》(http://slt.shanxi.gov.cn)《山西省生态环境公报》(https://sthjt.shanxi.gov.cn)《中国县(市)社会经济统计年鉴》^[18]及各市统计年鉴、省市县官网公报等，遥感影像及矢量数据来源于国家地球系统科学数据中心(www.geodata.cn)及中国科学院资源环境科学数据中心(http://www.resdc.cn)。

1.3 研究方法 1.3.1 预警指标体系构建

本研究遵循指标体系构建的科学性、可获取性及可操作性等原则的基础上，借鉴已有研究成果^[19-20]，从汾河流域本底特征出发，参考《山西省“十三五”环境保护规划》等指标要求，选取26项指标，通过DPSIR(驱动力—压力—状态—影响—响应)框架构建汾河流域土地生态安全预警指标体系(表 1)，其中，参评指标采用非零变换标准化统一处理^[21]，权重确定采用熵权法^[22]。

1.3.2 预警指数测算

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

${I_{ESW}} = \sum\limits_{i = 1}^n {{A_i}} \times {W_i} $

(1)

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

1.3.3 警情判别标准

结合汾河流域本底特征，参考已有研究设置本研究区的生态安全警情判别标准(表 2)，划定5级流域土地生态安全预警值^[20]。其中，土地生态安全预警值与研究区土地生态安全状况呈正相关关系。

1.3.4 空间自相关分析

空间自相关是用来分析研究区域范围内，各空间位置与各自周围邻近位置的同一属性值是否具有关联性，可以分为全局自相关和局部自相关^[23]。全局自相关用于揭示研究区整体某要素的空间相关性；局部空间自相关用于揭示研究区各子单元某要素空间相关性^[24]，为此，本研究采用LISA(local indicators of spatial association)聚类图分析土地生态安全警度在空间上的集聚特征^[25]。

$ {\rm{Moran's}}\;I = \frac{{\sum\limits_{i = 1j}^n {\sum\limits_{ij}^m {{W_{ij}}} } \left( {{x_i} - \bar x} \right)\left( {{x_j} - \bar x} \right)}}{{{S^2}\sum\limits_{i = 1j = 1}^n {\sum\limits_{ij}^m {{W_{ij}}} } }} $

(2)

式中：$ {S^2} = \frac{1}{n}\sum\limits_{i = 1}^n {{{\left( {{x_i} - \bar x} \right)}^2}} ,\bar x = \frac{1}{n}\sum\limits_{i = 1}^n {{x_i}} ;{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]。

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

$ \text { Moran's } I_{i}=\left(\frac{x_{i}-\bar{x}}{m}\right) \sum\limits_{j=1}^{n} W_{i j}\left(x_{i}-\bar{x}\right) $

(3)

式中：$ m = \frac{{\left( {\sum\limits_{j = 1,j \ne i}^j {x_i^2} } \right)}}{{(n - 1)}} - {\bar x^2}$ ；I_i值为“+”表示该空间单元周围相似值(高值/低值)的空间集聚；I_i值为“-”表示非相似值的空间集聚。

1.3.5 地理探测器

地理探测器是用来探测空间分异性及识别主导因子的空间统计学方法^[26]。由于汾河流域土地生态安全警度在空间分布上存在显著的异质性与分层性，且影响因子较多，因此，采用地理探测器识别主导影响因子，引入土地生态安全警度解释力q值来探测某因子X在多大程度解释属性Y(土地生态安全警度)，具体模型如下:

$ q = 1 - \frac{{\sum\limits_{i = 1}^L {{N_i}} \sigma _i^2}}{{N{\sigma ^2}}} $

(4)

式中：L代表因子X的分层(分类); N_i和N分别为层i和县区个数; σ_i²和σ²分别为层i和全局Y值的方差。解释力q介于[0, 1]；q值越大, 表明因子X对Y的解释力越强。

2 结果与分析 2.1 流域土地生态安全预警指数变化特征分析

2006—2018年汾河流域土地生态安全预警指数总体表现为先下降后增长的演变趋势(图 1)，2006—2012年，流域土地生态安全指数由0.383 4下降到0.337 0，平均下降幅度达12.10%；2013—2018年，流域土地生态安全指数下降趋势得到扭转，并以增长率3.18%逐渐提高，多年流域土地生态安全均处于较不安全状态，但2006—2018年拟合趋势表明流域土地生态安全由“较不安全”向“临界安全”状态转变。

从土地生态安全子系统层面对汾河流域土地生态安全预警指数变化特征来看，2006—2018年，驱动力与压力指数呈现平稳增长；2006—2012年状态指数增长了10.81%，2013—2018年下降了18.69%；影响指数多年基本维持在同一水平值，且出现小幅度阶段性波动；响应指数则呈现波动增长。分析发现，基于流域多年经济发展，特别是逐步提高的人均GDP及城镇居民可支配收入水平影响下，导致驱动力指数平稳增长。而经济发展有赖于土地系统高强度的开发与利用所提供的资源与环境，其中，化肥施用强度、土地垦殖系数及人均建设用地面积增长加剧了土地压力指数的提高。状态指数多年维持同一水平值限于城市人均公园绿地面积与第一、第二产业比重的正负效应驱动因子博弈的结果，也造成了影响指数与响应指数的关联性变化。

2.2 流域土地生态安全警度时空格局分析

由2006，2013，2018年预警指数时空分布发现(图 2)，基准期水平值介于[0.2, 0.6)，土地生态安全状态处于较不安全与临界安全，呈现中警和重警响应，土地生态安全整体水平偏低。

(1) 时间维度。2006—2018年流域土地生态安全水平经历“下降—提高”的过程，对比发现，除太原市区、汾阳市等19个评价单元土地生态安全水平呈向下开口“抛物线”的演变，其余县区土地生态安全水平均经历不同幅度的提高，演变趋势表现为向上开口“抛物线”增长。2006年，流域土地生态临界安全面积达16 043 km²，占流域总面积的39.47%，集中分布在流域中上游的太原、晋中盆地等县区，其中土地生态安全水平最高值0.53，最低值0.24，分别为太原市区与静乐县。2013年，流域“临界安全”同比下降至11.42%，处于临界安全范围的区县仅包括太原市区、寿阳县、汾阳县与古县，其他县区土地生态安全水平均为较不安全状态，最低值0.21出现在文水县。临界安全由2006年的11个县区下降至仅4个行政单元，到2018年，流域土地生态安全水平同比增长率为21.86%，达到临界安全水平的县区占研究区的23.80%，处于较不安全水平县区范围缩小11.42%，其中宁武县作为汾河源山区县增长率为33.55%，由较不安全水平提高到临界安全水平。综合来看，造成流域2006—2018年土地生态安全阈值的时间变化的是经济发展驱动下土地资源开发利用的反馈结果。流域社会经济的发展一方面为土地子系统提供了正向驱动力，如人均GDP与城镇居民人均可支配收入较2006年分别增长4.33%，43.01%，响应指标污水处理率等正向指标整体增长97.09%，为流域土地生态系统的良性循环提供了充足的经济技术等支撑。同时，高强度的土地开发与利用也加剧了流域土地生态安全水平的下降、警度的提高，另外，诱发土地生态安全水平下降的负向因子如化肥施用强度、土地垦殖系数、经济密度等对土地生态系统的施压存在滞后性，短期内表现的土地系统问题不显著，也由此带来治理响应的滞后，综合造成土地生态水平下降的累积结果。

(2) 空间格局。土地生态安全高值区集中分布于流域中上游东北部，零散分布于流域下游，整体呈现“大集聚、小分散”的格局特征。分析发现，经济发展程度、地形环境等是形成该格局的主要因子。经济较发达的太原辖区等，以省会城市政策倾斜、新型优势产业及环保投入等多项优势，使得土地生态水平均处于相对安全水平；而汾河中上游非完全山区县，自然环境较好，经济相对薄弱，人为对土地系统干扰力度较小，土地生态安全状态也能维持“临界安全”。土地生态安全的低值区主要位于汾河中下游盆地地区，人口经济高度集聚、产业结构不合理是形成该结果的主要原因。具体来看，第一类如尧都区等区县限于经济人口密度压力影响，产业以重化工业为主，土地生态安全状态“较不安全”；第二类如万荣县、稷山县等属于典型农业县，长期以来的高土地垦殖率及化肥使用强度等因素，使土地系统较长时间维持在“较不安全”状态。

2.3 流域土地生态安全警度空间关联特征分析 2.3.1 流域土地生态安全警度全局自相关

选用Moran’s I指数对汾河流域土地生态安全预警指数(ESWI)进行空间相关性分析(图 3)，2006，2013，2018年流域土地生态安全预警指数Moran’s I值分别为0.400，-0.070，0.373，总体经历“高—低—高”的下降—增长过程，表明流域土地生态安全警度经历“正相关—负相关—正相关”的转变过程。由此分析，汾河流域作为山西省人口经济高度集聚区，长期以来的高强度土地利用效率，远远超出其环境承载能力，累积的历史生态“欠债”多，土地生态安全形势严峻，环境问题多、治理难度大，周期长且波动显著，各评价单元差异性较大，造成流域土地生态安全警度全局空间相关性的变化。

2.3.2 流域土地生态安全警度局部自相关

局部Moran’s I指数与LISA聚类图可以分析局部单元与周边空间的关联性特征(图 4)，2006，2013，2018年汾河流域土地生态安全警度的局部自相关水平较低，主要显示为高—高(H-H)、高—低(H-L)、低—低(L-L)3类，2006年，H-H集中分布在流域中上游东北部太原市辖区、清徐县、晋中市辖区榆次区、寿阳县、太谷县5县，L-L分布在汾河源宁武县。H-H主要集聚在上述地区，是因为同处于太原都市圈影响范围，受经济发展程度影响，所提供的环保投资、三废处理率、绿化程度较高，另外，该年山西省实施《“十一五”环境保护规划》，全力推进蓝天碧水工程及“区域限批”，太原及晋中所辖区县取得较好的环境治理成果。而宁武县为L-L限于煤炭非法开采影响及汾河源早期旅游开发对环境保护评估不足，导致的环境破坏、土地生态的失衡。2013年，流域总体局部相关性水平出现下滑，H-H区域消失，流域多数评价单元相关性“不显著”，即评价单元与周围单元不存在相关。仅出现L-L在交城县、祁县，H-L在汾阳县，2006—2013年汾河流域总体社会经济发展加快，原有历史生态欠债并未还清且通过经济、人口密度加大等土地压力因子加重，使历史生态欠债进一步积欠，环境治理周期长，难度大。L-L分布的交城县、祁县是为地理条件及土地生产效率等制约形成局部低值单元；而H-L所在汾阳市为当年“绿色生态工程”实施的试点县之一，有效提高三废处理效率及绿化面积，也提高了土地生态安全水平。2018年，流域整体局部相关性提高显著，H-H区域扩大到汾河流域中上游8个区县，包括静乐县、阳曲县、古交市、太原辖区、清徐县、榆次区、太谷县及交城县。L-L分布在汾河流域下游襄汾县、翼城县、河津市，H-L分布在侯马市。2013—2018年汾河流域生态环境治理修复与保护工程取得阶段性成果，为汾河中上游流域(H-H集聚区域)土地生态安全系统改善起到了推动作用；而L-L区分布在汾河流域下游主要与长期以来基于地理条件的优越以此建立的高负荷土地利用模式所形成；H-L所在侯马市的产业结构转型升级以及环境阶段性治理成果对土地生态改善具有促进作用。

2.4 流域土地生态安全警度的影响因子分析

采用地理探测器识别流域土地生态安全警度主导影响因子(表 3)，使用GeoDetector前需对研究区进行格点化处理，划分5 km×5 km的1 626个格点，然后将所有格点所在位置的自变量及因变量提取作为输入数据^[27]；其次对输入数据进行重分类处理，将连续变量转化为类别变量^[28]。本文对26项影响因子X进行预处理，如地形坡度等依据生态环境状况评价技术规范(HJ192-2015)进行划分，其他影响因子均采用Jenks自然断点法分类，将影响因子X分为6类(层)。其中，q值为影响因子X对土地生态安全预警状况的解释力，q值介于[0, 1]，q值越接近于1，表示影响因子X对土地生态安全预警状况的解释力越强。由此得出，影响汾河流域土地生态安全的主导因子依次为X₃(人均GDP)、X₂₃(污水处理率)、X₂₅(建成区绿化率)、X₂₆(环保投资占GDP比例)、X₄(城镇居民可支配收入)、X₁₂(人均粮食产量)。

根据上述主导因子分析可知，汾河流域土地生态安全警度水平主要限于经济发展水平(X₃，X₄)、环境保护投入与治理力度(X₂₃，X₂₅，X₂₆)以及土地利用效率(X₁₂)等因素的影响。2006—2018年汾河流域经济发展水平稳步提高，人均GDP的增长、城镇居民可支配收入的提高，为流域土地科学利用提供了充足的资金支持，同时，也为土地生态安全状况的改善提供了正驱动力。环境保护投入与治理力度直接关系到土地生态安全子系统的循环与修复，2006—2018年，随着GDP的总体增长，流域内各县区也加大了对环境保护的投入与治理，对流域土地安全等级的优化作用显著。另外，土地的利用效率包括各土地利用类型单元内的增减更替，单位土地面积的GDP产值、粮食产量、能耗等。通过优化三生空间布局，产业结构转型升级以及科研技术投入等均可以提高土地利用效率。

3 结论

(1) 2006—2018年汾河流域土地生态安全预警指数随时间变化总体表现为向上开口的抛物线演变趋势，多年流域土地生态安全均处于较不安全状态，拟合趋势表明流域土地生态安全未来将由“较不安全”向“临界安全”状态转变。

(2) 流域土地生态安全警度时空格局，时间维度上，多年土地生态安全水平经历“下降—提高”的过程；空间格局上，土地生态安全高值区集中分布于流域中上游东北部，零散分布于流域下游，整体呈现“大集聚、小分散”的格局特征。

(3) 2006，2013，2018年流域土地生态安全警度的全局相关性经历“正相关—负相关—正相关”的转变过程。局部自相关主要显示为高—高(H-H)、高—低(H-L)、低—低(L-L)3类，高—高(H-H)主要分布在流域上中游，高—低(H-L)、低—低(L-L)主要分布在流域下游，与流域土地生态安全警度时空格局分布基本一致。

(4) 影响流域土地生态安全的主导因子为经济发展水平(X₃，X₄)、环境保护投入与治理力度(X₂₃，X₂₅，X₂₆)、土地利用效率(X₁₂)，所以，本文建议通过优化汾河流域生产、生活、生态空间布局，统筹全流域资源调配，注重流域上游源头保育、限制开发，中下游综合整治，建立流域生态补偿机制，加大环境保护与科研投入力度促进流域的生态安全水平的提高。