山区县土地生态脆弱性动态变化分析
——以青龙满族自治县为例

魏明欢, 胡波洋, 张贵军, 张蓬涛

(河北农业大学 国土资源学院, 河北 保定 071001)

生态环境是人类生存和发展的基础，但近几十年来，随着人口的快速增长以及工农业活动的不断加剧，生态环境受到越来越多的干预和破坏。日趋严重的生态环境问题使人类的生存和发展受到威胁，也引起了社会各界的广泛关注。2016年《国务院关于同意新增部分县(市、区、旗)纳入国家重点生态功能区的批复》(国函[2016]161号)中同意国家发展改革委提出的在国家重点生态功能区范围的基础上新增部分县(市、区、旗)纳入国家重点生态功能区，这表明国家对生态环境的保护力度进一步加大。

生态脆弱性研究是全球变化和可持续发展研究中的热点和重点[1]，到目前为止，生态脆弱性的研究已经形成了包括综合指数法、主成分分析法、景观评价法等在内的一系列分析方法，同时也形成了包括全国、区域、省市和流域在内的众多研究尺度[2]。封建民等[3]基于景观格局指数，从生态系统压力度、敏感度、恢复力3个方面构建指标体系，对榆阳区生态系统和各行政区的生态脆弱性进行分析和评价。廖炜等[4]基于RS和GIS技术，采用空间主成分分析法评价了丹江口库区的生态环境脆弱性，并进一步分析了土地利用变化与生态环境脆弱性的关系。张龙等[5]构建景观类型脆弱度指数计算模型和土地生态脆弱度指数计算模型，研究了黑龙江省宁安市1991—2010近20 a土地生态脆弱性的时空变化。总的来说土地生态系统脆弱性研究目前尚未形成一个统一的概念体系和分析框架，对山区县生态脆弱性时空分异特征以及演化规律的研究相对较少。

本研究以青龙满族自治县为研究区域，从自然条件和人为因素两方面选取7个指标构建评价指标体系，在综合指数法、变化斜率法等方法以及遥感、地理信息系统等技术手段的支持下，对青龙满族自治县2001年、2007年、2013年的生态脆弱性进行定量评价，揭示2001—2013年生态脆弱性的动态变化，以期为加强该县生态恢复及保护，促进资源环境与社会经济的协调发展提供参考，同时为山区县土地生态系统脆弱性研究提供借鉴。

1　研究区概况

青龙满族自治县地处河北省东北部边缘，地理位置为118.56°—119.64°E，40.08°—40.61°N。青龙满族自治县因其境内最大河流青龙河而得名，下辖25个乡镇，396个行政村，4个居民委员会，另有一个省级经济开发区；地形走势呈现东部、南部较低，西部、北部偏高，县境相对高度差1 766.3 m；气候属于大陆性季风气候，降水充沛；土壤类型以褐土、棕壤为主，两者约占总面积的97%；该县林果资源量多质优，全县林地面积4 920万hm2，森林覆盖率60%；受地貌特征影响，该县土地利用类型主要是以农业、林业、畜牧业和矿业用地为主。青龙满族自治县是地表水源涵养区，同时也是京津、冀东地区的生态屏障，因此该县生态与社会经济的协调关系着京津冀地区水资源、生态安全状况。但同时该县为山区县，昼夜温差大，降水多，土地冻融交替现象显著，水土流失、风蚀、水蚀现象普遍，生态环境极为脆弱，土地利用受到严重威胁。由于该县生态功能地位重要但土地生态系统极为脆弱，2016年，被国务院纳入国家重点生态功能区。

2　研究方法

2.1　评价指标体系的建立

参考已有研究成果[6-8]，结合青龙满族自治县实际，本研究从自然条件和人为因素两方面选取指标构建评价指标体系。

青龙满族自治县是典型的山区县，其生态环境脆弱性首先受到自然环境的影响。本研究从地形、土壤、气候、植被等方面选取指标表征自然环境对该县生态脆弱性的影响。坡度是水土流失、土壤侵蚀现象的重要影响因素，地形起伏度是表征地表单元地势起伏复杂程度的指标，两者都是地形因素的重要组成部分，因此本研究选取坡度和地形起伏度来表征地形对生态脆弱性的影响；植被是检验生态脆弱性的重要因素，NDVI因具有植被空间覆盖范围广、植物检测灵敏度高等众多优势条件，导致其在多种植被指数中应用最多也最广泛[9]，因此本研究选择NDVI表征植被对生态脆弱性的影响；土壤可蚀性可以定量地反映土壤对侵蚀作用的敏感性程度[10]，因此本研究选取土壤可蚀度表征土壤对生态脆弱性的影响；气候是生态环境形成和演化的重要动力，气候条件影响着生态环境的稳定性程度[11]，而降雨量的不均衡会产生引发水土流失的降雨侵蚀力，因此本研究选取年均降雨量表征气候对生态脆弱性的影响。人为因素对生态环境的影响是通过人类一系列的活动促进或减缓生态脆弱性的进程[12]。人口密度越大，说明人类活动越活跃，对自然环境的影响越强，因此人口密度是影响生态环境脆弱性的威胁因子；土地利用是人类活动作用于生态环境的直观表现[12]，不同的土地利用类型具有不同的生态功能，因此本研究选取土地利用类型表征土地利用状况对生态脆弱性影响。

2.2　数据来源与处理

2.2.1数据来源研究所需原始数据包括：(1) 遥感数据。本研究所需Landsat遥感影像数据来源于中国科学院计算机网络信息中心的地理空间数据云平台2001年、2007年的TM数据，2013年的OLI数据，成像时间为8月、9月份，分辨率均为30 m。(2) 数字高程(DEM)。来源于地理空间数据云平台分辨率为30 m的GDEM数据。(3) 土壤数据。来源于1∶30万青龙满族自治县土壤类型图以及《河北省土种志》。(4) 气象数据。来源于青龙满族自治县周边13个气象站点的实测数据。(5) 其他数据。包括2001年、2007年、2013年的青龙满族自治县国民经济统计数据以及青龙满族自治县行政区划矢量图。

2.2.2数据处理

(1) 土地利用类型数据获取。根据全国土地利用现状分类，并结合研究区实际情况，将青龙满族自治县土地利用方式划分为耕地、园地、林地、建设用地、水域、未利用地6类。利用ENVI 4.7软件对3期遥感影像进行预处理，然后采用监督分类和目视解译相结合的方法对遥感影像进行分类处理，并结合土地利用现状图，利用ENVI混淆矩阵工具进行精度检验，其中Kappa系数超过80%，解译精度达到研究要求[13]。

(2) NDVI数据获取。NDVI可以检测植被生长状况，计算公式如下：

(1)

式中：ρn为近红外波段的反射率；ρr为红外波段的反射率。

本研究利用ENVI 4.7软件的Transform-NDVI功能与ArcGIS掩膜提取工具，对遥感影像进行处理得到研究区3期NDVI数据。

(3) 坡度和地形起伏度获取。本研究基于青龙满族自治县DEM数据，利用GIS技术的Spatial Analyst Tools-Slope功能，获取青龙满族自治县坡度图，并利用移动窗口法获取其地形起伏度图。

(4) 气象数据的获取。利用研究区及其周围13个气象站点的实测数据，采用ArcGIS软件地统计分析中的反距离加权插值法进行插值得到研究区降雨数据。

(5) 土壤可蚀度数据获取。土壤可蚀度是评价土壤状况对土壤侵蚀中剥蚀、搬运等外应力敏感程度的指标，本研究根据青龙满族自治县土壤类型数据，并结合门明新等[10]采用EPIC模型对河北省部分土壤可蚀性的研究成果得到研究区土壤可蚀度分布图。

(6) 人口密度获取。根据青龙满族自治县2001年、2007年、2013年的国民经济统计资料获得研究区各乡镇三期人口数据，并利用ArcGIS软件赋值到青龙满族自治县行政区划图中，根据各乡镇面积计算得到人口密度，最后对行政区划图进行栅格化，分辨率为30 m，即可获得研究区三期人口密度栅格数据。

所有数据处理完毕后均为30 m×30 m的栅格数据。本研究以栅格为评价单元对青龙满族自治县三期生态脆弱性进行评价，并对其13 a间时空分异特征和演化趋势进行分析。

2.3　多重共线性检验

初选指标可能在一定程度上存在共线性问题，易造成评价结果的失真，因此在评价前应对评价指标进行共线性检验。本研究选取方差膨胀因子(VIF)作为共线性检验指标，VIF值越大，说明该变量与其他变量的相关性越大，通常选用VIF值等于5作为判断阈值[14]，当VIF值大于5时，说明所选指标的多元共线性较严重。

本文借助ArcGIS软件将青龙满族自治县划分为1 km×1 km的网格，并取网格中心点(3 706个)作为值提取点读取栅格数据，利用SPSS软件计算每个指标三期的VIF值(表1)，根据共线性检验结果可以看出，指标之间的VIF值均小于5，说明它们之间没有明显相关性。

表1　共线性检验结果

综上，本研究根据青龙满族自治县实际情况，同时兼顾指标综合性、独立性、代表性、可获取性等原则，从自然条件和人为因素两方面选取坡度、地形起伏度、土壤可蚀度、NDVI、年均降雨量、土地利用类型、人口密度构建生态脆弱性评价指标体系。

2.4　指标标准化及权重的确定

2.4.1指标数据标准化指标数据需要进行标准化处理，以消除变量的量纲。对于定量指标，本文采用极差标准化方法对原始指标数据进行标准化，公式为：

(2)

式中：Yij/是指第i个栅格第j个评价指标标准化值；Yij是指第i个栅格第j个评价指标原始值；Yimin是指第i个栅格第j个评价指标的最小值；Yiman是指第i个栅格第j个评价指标的最大值。具有正效应的指标包括坡度、地形起伏度、年均降水量、人口密度、土壤可蚀度，其值越大生态脆弱性越严重；具有负效应的指标包括NDVI，其值越小生态脆弱性越严重。

定性指标包括土地利用类型，本文参考相关研究[8,15]，通过赋值法对其进行定量化处理，之后利用极差标准化方法对其进行标准化，具体标准的选取情况见表2。

表2　分级赋值标准

2.4.2指标权重的确定熵权法是在综合考虑各因素信息量的基础上，根据各指标信息有序程度确定权重的一种方法，信息熵越小，指标权重就越大。作为客观综合定权法，熵权法能够消除人为因素的干扰，使评价结果更加科学合理[16-17]。

2.5　生态脆弱性评价

本文利用综合指数法，借助ArcGIS软件的Raster Calculator功能，将各指标标准化值与各指标权重的乘积之和作为生态脆弱性评价指数对青龙满族自治县生态脆弱性进行评价，其计算公式如下：

(3)

式中：EVIi表示第i个评价单元的生态脆弱性评价指数；Yij/表示第j项指标第i个评价单元的标准化值；ωj表示第j项指标权重。

在生态脆弱性评价指数的基础上，参考相关研究的评价标准[9,18]，并结合研究区自然地理特征，借助ArcGIS软件的Reclassify功能将青龙满族自治县生态脆弱性划分为5个等级，分别为微度脆弱、轻度脆弱、中度脆弱、高度脆弱和重度脆弱，划分标准见表3。

表3　生态脆弱性分级标准

2.6　生态脆弱性演化趋势分析

变化斜率法是利用最小二乘法对所有像元的变量与时间进行回归分析，预测其变化趋势，计算公式如下[19-20]：

(4)

式中：X为变化斜率；n为年数；EVIi为第i年生态脆弱性值。变化斜率为正，表明该区域生态脆弱性呈增加趋势，变化斜率为负，表明该区域生态脆弱性呈降低趋势。变化斜率显著性检验采用F检验，统计量计算如下：

(5)

3　结果与分析

3.1　研究区生态脆弱性评价

青龙满族自治县土地生态系统脆弱性评价指标权重表见表4，将权重结果带入式(3)中，得到研究区2001年、2007年、2013年生态脆弱性指数，进而得到3期土地生态脆弱性指数的空间分布图(图1)，根据等级划分标准，制作出生态脆弱性等级空间分布图(图2)，并在此基础上制作生态脆弱性等级分类表(表5)，以此分析研究区土地生态脆弱性时空分布特征。

表4　评价指标权重

图1　研究区生态脆弱性指数空间格局

3.2　研究区生态脆弱性时空分布特征

时间维度上，由图2和表5可以看出，2001—2013年微度脆弱区呈现先上升后下降的波动增加趋势；轻度脆弱区在研究期限内所占比重最大且呈现逐渐增加趋势；中度脆弱区呈现逐渐减少趋势；高度脆弱区呈现先下降后上升的波动减少趋势；重度脆弱区面积所占比重很小且在研究期限内变化不大。2013 年中度以上脆弱区面积为155 241.77 hm2，占研究区总面积的44.35%，说明研究区有将近一半区域处于土地生态脆弱性严重区。但从表5可以看出，2001—2013年青龙满族自治县中度以上脆弱性区所占比例由2001年的48.84%下降到2013年的44.35%，下降了4.49%，说明该县在13 a间积极进行土地生态系统的恢复和保护，并取得一定的成效，但由于山区自然条件的限制以及城镇化、工业化的发展，使得该县仍有部分地区生态条件较为恶劣。2001年、2007年、2013研究区生态脆弱性均以轻度脆弱性为主，面积为177 041.36 hm2，183 831.17 hm2，192 494.72 hm2，分别占研究区面积的50.56%，52.50%，54.97%，区域面积逐渐增大，中度以上脆弱区呈现逐渐减少趋势且大部分向轻度脆弱区发展，说明青龙满族自治县生态脆弱性朝着良好态势发展。

图2　研究区生态脆弱性等级空间格局

年份类型微度脆弱轻度脆弱中度脆弱高度脆弱重度脆弱2001年面积/hm22108．98177041．36163541．987442．760．06比重/%0．6050．5646．712．1302007年面积/hm25637．94183831．17155867．014797．711．31比重/%1．6152．5044．511．3802013年面积/hm22398．65192494．72149268．925972．850比重/%0．6854．9742．641．7102001—2007年面积/hm23528．966789．81-7674．97-2645．051．25比重/%1．011．94-2．20-0．7502007—2013年面积/hm2-3239．298663．55-6598．091175．14-1．31比重/%-0．932．47-1．870．3302001—2013年面积/hm2289．6715453．36-14273．06-1469．91-0．06比重/%0．084．41-4．07-0．420

空间维度上，研究区2001—2013年生态脆弱性普遍存在西北部高，中部、东部低的态势。高度脆弱区主要分布在西北部的青龙镇、八道河镇以及娄杖子镇，主要是因为该区域在县城周边，人口密度大、人类活动频繁，导致土地生态系统受到持续破坏。轻度脆弱区主要分布在研究区西南、中部的大部分区域，以及东南部的祖山林场，主要是由于该区域人为干扰程度低，地形较为平缓，同时该区域有包括桃林口水库、老岭自然保护区在内的众多生态环境保护重点区域，一定程度上减缓了生态脆弱性。

3.3　脆弱性演化趋势分析

依据变化斜率法以及显著性检验结果(图3，表6)，生态脆弱性显著提高的区域占全域面积的4.23%，且主要集中在轻度脆弱区，分布于茨榆山乡和大巫岚镇。可以看出该区域虽然有较低的脆弱性，但由于自然条件恶劣以及人类活动的影响，该区域土地生态系统的自我恢复能力逐渐减弱，生态环境的破坏程度逐渐加剧，因此应注意加强该区域生态环境的保护，防止生态脆弱性进一步恶化；脆弱性显著降低的区域占全域的10.06%，主要分布在八道河镇、龙王庙乡、三星口乡和干沟乡等，说明上述区域生态环境状况好转的趋势较为明显。由图3可以看出该区域多处于中度脆弱区和高度脆弱区，主要原因是随着水土保持、防护林建设等政策措施的实施，使区域土地生态系统逐步得到恢复；全域85.71%的像元不存在明显变化，且其中绝大部分处于中度脆弱区和轻度脆弱区，由此可以看出研究区大部分区域的生态脆弱性正处于中间过渡的重要阶段，如果该区域的生态环境保护能够得到足够的重视，及时采取修复和保护措施，生态脆弱性就会趋于减弱，生态环境状况就会好转；反之该区域的生态脆弱性会趋于增强，生态环境状况极易恶化。

表6　生态脆弱性变化趋势分析结果

图3　2001-2013年研究区生态脆弱性变化趋势及F值检验

4　讨论与结论

4.1　讨 论

(1) 本文所选指标重点突出了自然环境和人类活动对生态脆弱性的影响，指标选择具有普适性，数据具有可获取性。但所选指标能否较好地表征复杂的土地生态系统有待进一步验证。

(2) 本文借助遥感、地理信息系统等技术手段获取指标数据，但受数据来源和信息提取等方面的限制，数据的精度和质量有待进一步提高。

4.2　结 论

(1) 生态脆弱性时间维度分析表明研究区大部分地区处于生态脆弱性态势，但以轻度脆弱为主，且生态脆弱性朝着良好态势发展。

(2) 生态脆弱性空间维度分析表明：在研究期限内，研究区土地生态系统脆弱性普遍存在西北部高，中东部低的态势。高度脆弱区主要分布在西北部的青龙镇、八道河镇以及娄杖子镇，轻度脆弱区主要分布在研究区西南、中部的大部分区域，以及东南部的祖山林场。

(3) 生态脆弱性演化趋势分析结果表明：生态脆弱性显著提高的区域主要为轻度脆弱区，显著降低的区域主要为中、高度脆弱区，另全域绝大部分地区脆弱性变化趋势并不显著。

参考文献:

[1]徐广才,康慕谊,贺丽娜,等.生态脆弱性及其研究进展[J].生态学报,2009,29(5):2578-2588.

[2]张龙.宁安市土地利用/覆被变化及其生态脆弱性研究[D].哈尔滨:东北农业大学,2014.

[3]封建民,郭玲霞,李晓华.基于景观格局的榆阳区生态脆弱性评价简[J].水土保持研究,2016,23(6):179-184.

[4]廖炜,李璐,吴宜进.丹江口库区土地利用变化与生态环境脆弱性评价[J].自然资源学报,2011,26(11):1879-1889.

[5]张龙,宋戈,孟飞,等.宁安市土地生态脆弱性时空变化分析[J].水土保持研究,2014,21(2):133-137.

[6]韦晶,郭亚敏,孙林,等.三江源地区生态环境脆弱性评价[J].生态学杂志,2015,34(7):1968-1975.

[7]姚雄,余坤勇,刘健,等.南方水土流失严重区的生态脆弱性时空演变[J].应用生态学报,2016,27(3):735-745.

[8]马骏,李昌晓,魏虹,等.三峡库区生态脆弱性评价[J].生态学报,2015,35(21):7117-7129.

[9]刘世梁,田韫钰,尹艺洁,等.云南省植被NDVI时间变化特征及其对干旱的响应[J].生态学报,2016,36(15):4699-4707.

[10]门明新,赵同科,彭正萍,等.基于土壤粒径分布模型的河北省土壤可蚀性研究[J].中国农业科学,2004,37(11):1647-1653.

[11]黄方,刘湘南,张养贞. GIS支持下的吉林省西部生态环境脆弱态势评价研究[J].地理科学,2003,23(1):95-100.

[12]雷波,焦峰,王志杰,等.黄土丘陵区不同植被带典型小流域生态脆弱性评价[J].自然灾害学报,2013,22(5):149-159.

[13]许文宁,王鹏新,韩萍,等. Kappa系数在干旱预测模型精度评价中的应用:以关中平原的干旱预测为例[J].自然灾害学报,2011,20(6):81-86.

[14]孙瑜,史明昌,彭欢,等.基于maxENT模型的黑龙江大兴安岭森林雷击火火险预测[J].应用生态学报,2014,25(4):1100-1106.

[15]南颖,吉喆,冯恒栋,等.基于遥感和地理信息系统的图们江地区生态安全评价[J].生态学报,2013,33(15):4790-4798.

[16]熊黑钢,付金花,王凯龙.基于熵权法的新疆奇台绿洲水资源承载力评价研究[J].中国生态农业学报,2012,20(10):124-129.

[17]赵丽,朱永明,付梅臣,等.主成分分析法和熵值法在农村居民点集约利用评价中的比较[J].农业工程学报,2012,28(7):235-242.

[18]徐涵秋.水土流失区生态变化的遥感评估[J].农业工程学报,2013,29(7):91-97.

[19]陈操操,谢高地,甄霖,等.泾河流域植被覆盖动态变化特征及其与降雨的关系[J].生态学报,2008,28(3):925-938.

[20]陈佳,杨新军,尹莎,等.基于VSD框架的半干旱地区社会—生态系统脆弱性演化与模拟[J].地理学报,2016,71(7):1172-1188.