基于自然-经济-景观多源体系的吕梁山土地生态风险评价

宋晓静，周淑琴，王 勇，荆耀栋，刘立文，王贝贝

（山西农业大学资源环境学院，山西 太谷 030801）

生态环境评价是一种能有效支持生态系统管理的新模式［1］，土地生态风险评价与生态系统服务价值的测算是当前生态环境评价的组成部分［2，3］。土地生态风险评价以景观生态学的生态过程与空间格局的耦合为视角，注重风险的时空异质性和尺度效应，以实现多源风险受体的综合表征及其空间可视化［4］。当前，构建景观干扰度和脆弱度指数［5，6］揭示区域生态风险分布特征，已成为评价区域生态环境质量的新方法之一。刘珍环等［7］通过脆弱度和干扰度指数评价三江平原的景观生态风险，揭示该区域生态风险的空间分布与格局变化。生态系统服务理念基于自然资本及其生态服务功能［8］，以期实现人类福祉最大化［9］，包括供给、调节、支持和文化服务4方面［10-13］。由于自然资本在城镇扩张中不断遭受破坏和侵蚀，导致自然景观转为人为景观，生态环境调节功能衰退，“碳汇”功能减弱或丧失［14］，生态系统服务也在不断下降［15］。马倩等［16］采用生态系统服务价值计算公式与单位面积价值当量表，分析新疆艾比湖生态脆弱区1990—2007 年和生态系统服务价值的变化对土地利用特征的影响。

为实现多视角下的生态环境评价，现阶段学者们的研究从单一对生态风险评价或生态系统服务价值测算逐步走向二者的融合［17］。康鹏等［18］将生态系统服务引入生态风险评价，建议结合生态系统的特征及评价要求，建立动态评价模型，模拟不同情景下生态风险高低对生态系统服务的影响，有效改进传统评价模型。曹祺文等［19］以生态系统服务退化为基础构建多维度指标体系，建立基于生态系统服务的景观生态风险评价（ESRISK）框架，创新了生态风险评价方法。钟小龙［20］分别计算景观指数和土地生态系统服务价值，并对二者关系进行相关性分析，结果证明，二者之间存在协同作用，景观指数的变化对土地生态系统服务价值有一定的影响。

十九大以来，习近平总书记提出“加强生态文明体制改革，建设美丽中国”。生态文明建设成为热点，面向生命共同体治理的生态环境评价显得尤为重要。吕梁山作为拱卫京津冀和黄河生态安全的重要屏障［21］，区域内山、水、林、田、湖、草、沙各自然要素构成生命共同体。区域生态环境复合评价体系包括自然、经济和景观生态安全［9，22］，前期研究多选取某层或二者结合展开。鉴于此，本研究以人为活动频繁的自然地理单元——吕梁山为例，以“地理+生态”［23］的模式剖析土地利用变化下的生态效应，从植被退化、生态系统服务价值、景观生态角度3 个层面展开评价，深度挖掘吕梁山区生态风险及问题，提出系统保护与修复建议。这是加强生态文明建设的现实需求［23］，也是实现区域山水林田湖草沙生命共同体的治理修复，提高该地区生态风险预警和管控能力的重要参考。

1 研究区概况

吕梁山位于山西省西部，黄河山西段东岸。吕梁山以西为水土流失和风沙区，以东为盆地与山区的过渡带，山脉主脊两侧多分布为天然林地（图1）。海拔高度在543～2 772 m，大部分区域集中于1 100～1 800 m。以山地和丘陵为主，地形起伏大、地貌类型多样，林地、草地在农用地中占有重要位置，坡耕地多且布局零散。大陆性季风气候，降水分布不均且多集中于6—9 月。土壤类型以黄绵土、褐土、浅栗褐土为主，著名景点关帝山、芦芽山等分布有成片原始林地，森林资源丰富，华北落叶松、油松、云杉、白杄和青杄林林地呈条状自北向南分布。西部集中分布有成片发育良好的亚高山草原，为省内重要牧场。吕梁山作为国家著名革命老区，是拱卫京津冀和黄河生态安全的重要屏障，也是山西省生态建设和脱贫攻坚的主战场，在发展过程中坚持“以创新促增绿，以改革促增收”的目标，力争改善老区人民生产生活条件，打好打赢生态治理和脱贫攻坚两场战役。

图1 研究区位置

2 数据源与研究方法

2.1 数据源及预处理

数据来源于中国科学院资源环境与数据中心（http：//www.resdc.cn）2000 年、2010 年、2018 年土地利用动态监测数据，并按照国家土地利用现状分类标准（GB/T21010—2017），结合山西省实际情况，将土地利用类型分为耕地、林地、水域、草地、建设用地和未利用地6 类。高程及行政区划矢量数据下载自地理空间数据云（http：//www.gscloud.cn）。NDVI 数据来源于MOD13Q2 16d 合成产品，空间分辨率500 m，时间序列2000—2018 年，采用像元二分模型合成研究区逐年植被覆盖度（FVC）数据。通过重采样将数据统一成30 m。所涉及到的矢量行政区数据来自于地理国情监测云平台（http：//www.dsac.cn/）。其他基础数据来自历年中国统计年鉴、山西省统计年鉴汇编，相关数据经过基础数据计算得到。

2.2 研究方法

2.2.1 土地利用生态风险指数核算模型 根据研究区范围和平均板块面积（2～5 倍），采用5 km×5 km 格网进行网格化采样，将吕梁山划分为1 840 个生态风险小区，计算每一小区内的生态风险指数值（ERI），并将其作为样本中心点的生态风险水平。引入景观干扰度Ui、脆弱度Li和损失度Ri来计算景观生态风险指数［4，5］。

式中，ERI为研究区景观生态风险指数，Ak是景观总面积，Aki是某一景观类型i的面积，m为景观类型的数量，Ri是第i类景观损失度指数。

景观干扰度Ui由景观破碎度Ci、景观分离度Fi、景观优势度Di赋予权重后计算得出，结合前人研究成果［4，5，22］，将其权重分别赋予 0.5、0.2、0.3。景观脆弱度用来表示不同景观类型遭受到外界干扰后的敏感程度，根据吕梁山土地利用分类情况并结合已有成果，将景观脆弱性分为6 级，从高到低依次为耕地、草地、水域、未利用地、林地、建设用地，进行归一化处理。

2.2.2 生态系统服务价值核算模型 基于当量因子表以及谢高地等［23-25］在2015 年修订和补充的“单位面积生态系统服务价值（ESV）当量表”，定义为每公顷农田生态系统中粮食生产净利润作为标准下ESV的当量。吕梁山地区的农业用地中，耕地以旱地为主；林地以有林地、灌木林地和针阔叶混交林为主。本研究尺度均在山西省内，区域内部自然地理特征基本一致，农业生产都以旱作农业为主，各生态系统属性相近，可以利用传统方法核算生态系统服务价值（表1）。

表1 2000—2018 年吕梁山生态系统生态服务价值均值 （单位：108元）

式中，ESV、ESVf分别为研究区第i类服务各功能的价值和生态系统服务总价值；Ai为第i种土地利用类型的面积；VCi为第i种土地利用类型的生态系统服务价值系数。

2.2.3 植被退化趋势评价 根据GB19377—2003《天然草地退化、沙化、盐渍化的分级指标》［26］、吕梁山植被退化现状和相关研究，选择植被覆盖度（FVC）作为植被退化监测指标，将吕梁山2000—2018 年植被退化程度分为3 级：未退化、轻度退化、中重度退化（表2）。采用Sen+Mann-Kendall 趋势检验法对2000—2018 年FVC退化趋势以及空间格局进行分析。Sen+Mann-Kendall 作为一种非参数统计检验方法，样本无需遵从一定分布，规避异常值干扰的能力强，广泛用于检验时间序列变化趋势的研究［27］。显著性检验公式如下。

表2 基于FVC 的植被退化程度分级

式中，n代表时间序列长度，S代表时间序列中对偶值的大小关系，Var(S)代表S的方差，Zc为检验结果。一般在显著性水平α=0.05 下进行统计检验，|Z|≥ 1.96 呈显著变化。

趋势度计算公式如下。

式 中 ，FVCa、FVCb代 表 第a年 和 第b年 各 像 元FVC的最大合成值。β值为正，表示FVC呈上升趋势，负值则表示为下降趋势。该部分通过R编程实现。

2.2.4 双变量空间自相关模型 双变量空间自相关模型是探索性空间数据分析（Exploratory spatial data analysis，ESDA）的一种。ESDA 是通过描述某一地理事物的空间分布，并对其进行可视化处理，识别其空间布局、集聚特点和异常，揭示对象的空间模式和作用机制［28］。ESDA 包括全局空间自相关和局部空间自相关。其中，以全局空间自相关研究整体空间关系，以局部空间自相关研究空间的变异状况。通过ArcGIS 和GeoDa 空间计量分析软件引入ESDA，构建双变量探究生态风险指数与生态系统服务价值的相互作用机制。

1）全局空间自相关。全局空间自相关通常采用Moran’s I 指数表示。Moran’s I 的绝对值越接近 1，表明研究单元的空间自相关程度越强，计算公式如下。

式中，I为生态风险指数i与生态系统服务价值j的双变量自相关系数，yi、yj表示某个单元的生态风险指数和生态系统服务价值，n为研究单元的数量，代表两个指数的均值；ωij是相邻权重，σs、σr代表两个指数的方差。

2）局部空间自相关。LISA 图用以识别生态风险区高值、低值的空间关联模式，计算公式如下。

Moran’s I 值一般在［-1，1］区间内，小于 0 为负相关，大于0 为正相关，等于0 则表示不相关，S2为样本方差。根据LISA 图获取空间分布关系，并依据其空间关联性划分5 种集聚类型：高高聚集（H-H）、高低聚集（H-L）、低高聚集（L-H）、低低聚集（L-L）以及不显著。

3 结果与分析

3.1 植被退化趋势评价

受气候变化、人类活动以及国家生态保护措施的多重影响，2000—2018 年吕梁山植被覆盖度呈逐渐增加趋势。趋势系数介于（-0.034 7～0.059 4）/19年，呈下降和上升趋势的区域面积分别占研究区总面积的21.90%和78.10%，显著相关（P＜0.05）的区域面积占总面积的84.94%，吕梁山植被变化多具有显著性，其中显著不增加区域面积占比32.35%，显著增加区域面积占比40.70%，不显著不增加与不显著增加区域面积占比总和为0.065%（图2）。从整体来看，植被改善以及基本稳定不变的面积高达73.115%，区域植被治理与保护的优势有所显现。但也有部分地区出现植被中重度退化现象，集中在吕梁山西部的临县、兴县、柳林、石楼和北部河曲县、偏关县及中部娄烦县以及静乐县等地。这些地区的共同点之一是长期处于贫困与半贫困状态，为谋求发展所需的生态系统服务数量和种类越来越多以及产业优势的缺乏导致这些县（区）长期以牺牲林草地资源为代价发展经济，对自然资源的利用已超越生态系统本身提供限度［9］，尤其以近几年脱贫任务艰巨，经济发展对生态环境的破坏更为突出。

图2 植被退化趋势与相关性检验

3.2 生态系统服务价值时空变化分析

依据表3，吕梁山2000—2010 年生态系统服务价值有所增长，总价值从507.803×108元增长至507.874×108元，共增长 0.014%；2010—2018 年生态系统服务价值出现递减的趋势，总价值从507.874×108元减少至 499.761×108元，共减少 1.597%。不同用地类型的生态服务功能和价值不同，一般来说建设用地生态服务价值低，而农用地、水域的生态服务价值高。研究区各土地利用类型的生态系统服务价值占比大小依次为林地、草地、耕地、水域、未利用地、建设用地。水域和林地、草地服务价值增加最显著。林地、草地、耕地是吕梁山生态系统服务价值的主要贡献类型，20 年间林地均保持在50%以上的贡献率，而耕地的贡献率逐年递减；草地的贡献率则呈现先增加后减少的趋势。近10 年来，工矿开采、天然林木砍伐以及公路、铁路建设等人类活动剧烈，对当地生态环境系统造成严重干扰，也诱发了许多自然灾害和环境问题。建设用地扩增的同时，农业用地的面积被逐渐蚕食［29］，与之配套的各类生态系统服务价值逐渐减小，建设用地带来的生态系统服务负面价值［30］越来越明显。

表3 吕梁山2000—2018 年生态系统服务价值

结合ArcGIS10.2 Natural Breaks 和研究区实际状况，对生态系统服务价值进行可视化表达，并将其划分为 5 级：低［0，900）、较低［900，1 800）、中［1 800，2 900）、较高［2 900，4 100）、高［4 100，+∞）（图3）。吕梁山单位面积生态系统服务价值以较低、中等和高等3 个等级为主，面积平均占比为30.43%、21.58%和18.91%；高值区面积基本保持不变，管涔山北段、南段、云中山、中部吕梁山最为稳定，其他地区则零散分布。2000—2010 年低值区显著减少，较高、中等、高值区的面积比例则有所提高，改善优于退化，生态服务价值有所提高。2010 年之后，除低值区、较低值区域面积有增加外，中等、较高值区面积均有所减少，导致整个区域的生态系统服务价值有所下降，而吕梁山脉北段、中段以及南部火焰山因其特殊的地理区位，森林资源所具有的水源涵养、气候调节等优势［31，32］，生态系统服务价值最高，变化也最稳定。

图3 2000—2018 年吕梁山各项生态系统服务价值

3.3 景观生态风险时空变化分析

为分析吕梁山景观生态风险的空间分异特征，利用ArcGIS10.2 的Natural Breaks 将研究区的景观生态风险划分为5 级，赋值区间为高风险区（ERI＞1.00）、较高风险区（0.16＜ERI≤1.00）、中风险区（0.12＜ERI≤0.16）、较低风险区（0.09＜ERI≤0.12）、低风险区（ERI≤0.09）。2000—2018 年吕梁山景观生态风险整体呈现上升趋势，其中较低、中生态风险区占比较大，高风险区生态风险占比较小（图4、图5）。

图4 2000—2018 年吕梁山景观生态风险各等级面积占比

图5 2000—2018 年吕梁山景观生态风险空间分布

2000 年吕梁山低、较低生态风险分布区域分别占研究区总面积的21.30%、23.64%。低风险区集中分布于中部山脉周边，以中阳县为主，北部管涔山周边以兴县、岢岚县为主，南部火焰山附近以吉县、乡宁县和蒲县为主。较低风险区主要分布在低风险区周边，以云中山附近居多。这两类风险区的土地利用类型以林地为主，用地类型单一且集中连片，由于林地的景观脆弱度较低，该区域生态风险等级较低。中风险区主要分布在耕地、草地等景观相间分布区，占总面积的52.45%。较高、高风险区主要分布在娄烦、古交等靠近省会城市的区域，分别占总面积的2.45%、0.16%。该区域有耕地分布，以坡耕地为主，人类活动强度的增加造成景观破碎化、林草退化、农用地转为建设用地等生态环境问题。

2010 年吕梁山生态风险的空间分布规律较2000 年并未出现明显变化。其中低、较低等级所占面积较2000 年增加了1.36%，面积变化的区域主要是低、较低等级继续向外扩展，面积增大。这些地区退耕还林保护强度增大，大面积农业景观得以保存。中等级分布面积小幅度减少，较高、高风险区的面积没有明显变化。说明这10 年间，吕梁山所在辖区的市（县）统筹推进荒山造林、水源保护、退耕还林和草食畜牧，生态状况保持稳定。

2018 年吕梁山较高、高风险区的面积及空间分布均出现显著变化，较2010 年分别增加了5.54%、0.05%；低、较低风险区较2010 年相比减少2.55%，中风险区与前期相比减少3.04%。较高、高生态风险区除了分布在靠近省会城市的县（区）外，建设集中化使高风险区有向外蔓延的趋势，造成区域内生态保护压力增大。另外，部分县刚完成脱贫任务，经济产业转型升级的需求导致大量农业类景观向建设用地转换，景观破碎度、脆弱度也在持续扩大。

3.4 生态风险与生态系统服务价值空间相关性分析

将生态风险指数定义为变量1，生态系统服务价值定义为变量2，定义空间权重后，带入Geoda 中进行双变量空间自相关分析。结果表明，生态风险指数与生态系统服务价值成负相关，即生态风险程度越高，相对应的生态系统服务价值越低。3 期的Moran’s I指数分别为-0.179、-0.151、-0.150，全局自相关指数逐年提高，聚集趋势有所增强。

研究区生态风险指数与生态系统服务价值的局部空间自相关LISA 分析结果如图6 所示，共划分为5 种模式，不显著、高风险-高价值、低风险-低价值、高风险-低价值及低风险-高价值。结果表明，相关性显著区域代表着ESV和ERI呈现出一定的空间集聚状态。主要分布在北部管涔山、云中山周围以及中部吕梁山、南部火焰山的大部分以及山区边缘区域的县（区），包括北部的偏关、保德以及中部吕梁市的兴县、临县、柳林，南部临汾市的汾西、洪洞等县。不显著区域代表着ESV和ERI空间集聚不明显。该部分面积占研究区总面积的70%左右，以零散分布的耕地及草地为主，ERI和ESV也都呈现均匀分布，未出现明显的异常值。高风险-高价值代表研究区ERI较高、ESV较高的区域。集中在吕梁山脉周边的小部分地区，方山、交城县周边集中分布有大片优质草原，包括西华镇草原、云顶山草原及庞泉沟国家级自然保护区等，地区在不断开发旅游业扩宽产业结构的同时，也带来一定程度的生态风险。但绝大部分区域仍处于低风险-高价值，即ERI较低、ESV较高的区域，这也与实际情况相符。吕梁山以林地景观为主，植被覆盖率高，景观优势度较高且分布集中，支持、供给、调节服务均能得到较好发挥，较少受到人为活动干扰，生态风险指数也偏低。高风险-低价值代表研究区ERI较高、ESV较低的区域。该区域的形成与当地的资源禀赋有关，西部的保德县、兴县、临县等都是国家级贫困县，产业结构的失衡与当地较低的经济水平导致当地的生态系统服务价值偏低，生态风险偏高。低风险-低价值代表研究区ERI较低、ESV较低的区域。该区域分布在人为干扰较强的县（区），所提供的生态系统服务价值以满足人类需求为主，而人为干扰对环境的不利影响使得该区处于相对稳定的状态，风险指数也相对较低。

图6 2000-2018 年吕梁山景观生态风险局部自相关

研究从自然、经济、景观3 个维度对吕梁山地区的生态环境状况展开评价。结果表明，以北段管涔山、云中山，中段吕梁山，南段火焰山为代表的主体山脉周围，以林地景观为主，整体植被覆盖度最高，生态服务价值的高值、较高值区也多集中分布于此处，受人类活动干预较少，生态环境稳定，良好的生态环境带来区域内较低的生态风险。而山脉周围的地区，尤其以省会、地级市集中连片区以及刚完成脱贫任务的部分县（区）而言，受中心城市辐射以及产业转型的影响，建设用地扩张侵占了大量的农业用地，带来地区植被退化、生态系统服务价值的降低以及生态风险区的集聚，如果人类继续不合理地开发利用资源，逆向演化还会持续扩张。区域内实施的林草地封育保护、水源地保护、防风固沙等措施，使得前10 年间局部生态环境有所改善，生态风险降低，但由于生态保护效应的滞后性［33］及生态环境的脆弱性，山区的生态功能短时期内难以恢复，总体仍在恶化，生态风险呈增高趋势。因此，应当立足全局，以吕梁山为中心，妥善处理当地人口与水资源、土地资源和生态环境之间的关系，制定科学合理的规划，以实现生态系统的良性循环和可持续发展。

4 结论与建议

4.1 结论

2000—2018 年吕梁山植被覆盖度呈增加趋势。植被改善区以及基本稳定区的面积占比73%以上，区域植被治理与保护的优势有所显现，以贫困地区为代表的局部区域出现明显的植被退化趋势。

19 年间吕梁山生态系统服务价值先增后降。农用地是吕梁山生态系统服务价值的主要贡献类型；从各生态系统服务指标来看，土壤形成与保护和生物多样性保护的生态价值最高，食物生产与娱乐文化的生态价值占比较低。研究区生态系统服务价值以较低、中等和高等3 个等级为主，以山脉为代表的高值区变化较稳定。

研究区景观生态风险整体呈上升趋势，较高、高风险区的基质景观以耕地为主，受人为干扰后损失度变高。受人为因素的影响，土地利用类型的转化强度也在持续扩大，景观类型转换区的脆弱度和破碎度都有所增强。

生态风险指数与生态系统服务价值成负相关，Moran’s I 指数逐年提高，聚集趋势有所增强。5 种分类模式中，不显著区域面积占总面积2/3 以上，高风险-高价值、低风险-高价值集中分布在主体山脉及其边缘，高风险-低价值地区以常年贫困地区的资源禀赋为代表，低风险-低价值地区面积最小，处于相对稳定的状态。

4.2 建议

吕梁山景观生态安全变化中，“源景观”［34］（主要包括自然保护区、森林公园等地）和景观廊道（河流和人工湖等）在城市化及资源整改转型中受到的胁迫相对较小，保持较为稳定、和谐的景观结构和生态服务功能，生态风险指数的变化是由于耕地等农业景观逐步减少，建设用地不断增加以及工业用地的扩张造成的。为保持区域生态可持续发展，需要持续保护生态源地和景观廊道，加强对国家级自然保护区的管护［35］。因地制宜，合理利用耕地，修建梯田，提高植被覆盖率，多措并举治理水土流失，有效降低人类活动对自然景观的威胁；同时，以新一轮土地利用总体规划和功能区划为依据，增加林地保护力度，严格控制林木采伐量和林地占有量［36］，加强水源涵养林建设，增加生态用地数量，改善生态服务功能，保护生物多样性。考虑到建设用地开发对景观生态安全的影响，未来研究区规划中还需进一步优化建设用地的空间布局，分区域分行业控制建设用地密度，增加建设用地的绿化空间，有效缓解生态环境的压力；通过调整建设用地空间架构，合理规划城乡建设用地，有效应对受人类与自然过程扰动的景观破碎化地区的严峻环境和发展问题，在保证生态系统稳定性的同时稳步提升生态系统服务价值，以提高全域生态安全。