结合主导生态功能与生态退化程度的生态修复分区方法

摘要

本发明涉及一种结合主导生态功能与生态退化程度的生态修复分区方法，获取待研究区的原始数据，对原始数据进行预处理和信息提取，绘制底层数据图；计算水源涵养、防风固沙、土壤保持和生物多样性维持这四个生态功能要素；并计算生态退化程度，生成生态退化程度图；分三级划分生态修复分区，生成各级分区图；将各级分区图进行叠加，生成最终的生态修复分区图。本发明所述方法克服了以往生态功能区划只重视生态系统主要服务功能空间布局的局限性，实现了生态修复分区的定量划分，提高了生态修复工程的针对性和精确性，节省了人力、物力。

说明书

技术领域

本发明涉及生态恢复领域，尤其涉及一种结合主导生态功能与生态退化程度的生态修复分区方法。

背景技术

随着经济社会发展，资源开发利用与生态环境保护矛盾突出，表现为生态环境敏感性越来越大，而人类对生态系统服务功能的需求越来越多，因此合理安排生态恢复和经济发展布局，科学高效的开展生态修复成为关注焦点。生态修复分区理论基于生态区划，而生态区划的概念最早是由俄国地理科学家道库恰也夫于1899年提出的，1905年英国生态学家赫伯特森(Herbertson)针对全球尺度上的自然区域单元开展了区划研究，成为生态区划研究的重要里程碑。此后，生态区划受到越来越多的科学家关注，但早期的区划只是围绕自然生态系统植被展开，主要考虑水热因子对生态系统的影响，建立了气候指标体系，用以区分不同类型的生态系统。随着研究的进一步深入，各国生态学家在基于自然因子的生态区划方面取得了显著的进展，1976年美国生态学家贝利(Bailey)第一次提出了较为完整的生态区划方案。随着经济社会的发展，人口急速膨胀，出现一系列生态环境问题，区划方法也将人为因素的影响纳入到了区划指标中。80年代以来，我国各单项区划和综合区划日趋完善，从生态胁迫、生态敏感性、生态压力角度构建进行地理空间分区，并从不同省域、城市、县域、流域和湖泊等不同区域尺度筛选指标进行生态功能区划，结合经济社会发展趋势分析区域内生态环境问题及其形成机制。为了进一步保护生态系统，2008和2011年环境保护部发布了《全国生态功能区划》和《全国主体功能区规划》，在生态环境现状、生态环境敏感性、生态功能重要性评价的基础上划定全国生态功能区，对维护区域生态安全具有重要意义。

生态修复分区是开展生态修复的前提，目前生态分区主要集中在生态功能区划，且生态修复分区存在主观性强，难以量化、不易操作的缺点，从而使生态修复工程针对性较差，造成生态修复中的资源浪费，无法合理高效地利用生态修复资源来提升当地生态功能，因此需要一种系统的可量化的生态修复分区技术方法。

发明内容

鉴于上述的分析，本发明旨在提供一种结合主导生态功能与生态退化程度的生态修复分区方法，用以解决现有分区方法主观性强，难以量化、不易操作的问题。

本发明的目的主要是通过以下技术方案实现的，提供一种结合主导生态功能与生态退化程度的生态修复分区方法，包括：

步骤S1.获取待研究区的原始数据，对原始数据进行预处理和信息提取，绘制底层数据图；

步骤S2.计算生态功能要素，并生成各生态功能要素栅格图，所述生态功能要素包括水源涵养、防风固沙、土壤保持和生物多样性维持；

步骤S3.计算待研究区的生态退化程度，生成生态退化程度图；

步骤S4.分级划分生态修复分区，生成各级分区图；

步骤S5.将各级分区图进行叠加，生成生态修复分区图。

其中，原始数据包括行政区划数据、高程数据、气象要素数据、土地覆盖类型数据、遥感影像数据。

对原始数据进行预处理和信息提取，绘制底层数据图进一步包括：

对行政区划图数据进行地图配准，再进行边界矢量化，得到待研究区的矢量边界图；

对高程数据进行配准，然后裁剪，再提取待研究区出的坡度、坡向以及坡度坡长因子，并绘制坡度图、坡向图、坡度坡长因子图；

对气象数据，包括气温、降水、风速、风向，分别进行插值，生成各气象要素的栅格数据，再进行裁剪；

对土地覆盖类型数据进行配准，再进行裁剪，获取待研究区的土地覆盖类型图；

针对遥感图像，选取红外波段与近红外波段计算植被归一化指数NDVI，在此基础上，运用像元分解模型，获取植被覆盖度数据，经过裁剪后，绘制待研究区的植被覆盖度图。

步骤S2中，水源涵养量的计算采用降水贮存量法，生成表征水源涵养量的图；防风固沙是基于植被覆盖度图，绘制防风固沙量的空间分布图，生成表征防风固沙量的图；土壤保持采用通用水土流失方程进行计算，生成表征土壤保持功能的图；生物多样性维持采用Invest-biodiversity模型进行计算，模型中采用生境质量指数进行评价，最终生成表征生物多样多样性维持功能的图。

步骤S3中，生态退化程度是依据理论NDVI与实际NDVI比值，并参照退化等级划分标准，确定生态退化等级，生成生态退化等级图。

步骤S4中，分级划分生态修复分区进一步包括：根据立地条件划分生态修复分区的一级区，生成生态修复一级区图，所述立地条件包括气候气象、地形地貌、植被覆被；根据主导生态功能要素划分二级区，生成生态修复二级区图；根据生态退化程度划分三级区，生成生态修复三级区图；

其中，所述主导生态功能要素的确定方式为：

将水源涵养、防风固沙、土壤保持和生物多样性维持这四个生态功能要素的计算值进行标准化，再将标准化后的值进行等级划分；在各个栅格位置上，将上述四个生态功能要素标准化值中，最大等级对应的生态功能要素作为主导生态功能要素，并根据各栅格确定的主导生态功能要素进行对应赋值。

上述根据生态退化程度划分三级区，进一步包括：按照生态退化程度，将待研究区划定为未退化、轻度退化、中度退化、重度退化这4种退化等级；将轻度退化划定为预防修复区，中度退化划定为调控修复区，重度退化划定为重点修复区；其中，在未退化区进一步采用自然断裂法，按照实际NDVI值从大到小分为3个等级，所述3个等级是优先保护区、重点保护区和一般保护区。

本发明有益效果如下：构建生态修复分区的生态功能要素，确定主导生态功能要素；并结合生态退化程度判定生态修复难易程度；分3级划分生态修复分区，实现了定量划分，克服了以往生态功能区划只重视生态系统主要服务功能空间布局的局限性，提高了生态修复工程的针对性和精确性，节省人力、物力。

本发明的其他特征和优点将在随后的说明书中阐述，并且，部分的从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。

附图说明

附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。

图1为生态威胁因子属性；

图2为土地覆盖类型对威胁因子的敏感度；

图3为划定各级生态修复空间分区的流程图；

图4为生态修复分区的要素；

图5为实施例的生态修复分区结果。

具体实施方式

下面结合附图来具体描述本发明的优选实施例，其中，附图构成本申请一部分，并与本发明的实施例一起用于阐释本发明的原理。

本发明的实施例以内蒙古巴林右旗作为生态修复分区的案例，分级进行了生态修复分区划定，具体步骤如下。生态修复分区中的数据处理均通过地理信息系统软件进行，本实施例具体使用但不局限于ArcGIS软件。

1.获取巴林右旗研究区的行政区划数据、高程数据、气象要素数据、土地覆盖类型数据、遥感影像等原始数据，并进行预处理和信息提取，并绘制巴林右旗研究区的各个底层数据图。

1.1.获取巴林右旗行政区划图，将该行政区划图作为底图，首先进行地图配准，再进行研究区域的边界矢量化，得到巴林右旗矢量边界图。该步骤可以使用ArcGIS软件中的Georeferncing模块进行地图配准。

1.2.获取覆盖巴林右旗区域的DEM(Digital Elevation Model)数据，首先进行配准，然后利用巴林右旗矢量边界，裁剪DEM数据，得到以上述边界界定的DEM；在经过裁剪的DEM数据基础上，进一步提取出巴林右旗研究区的坡度、坡向以及坡度坡长因子，并绘制坡度图、坡向图、坡度坡长因子图。该步骤可以在ArcGIS软件中的surface analysis模块中进行。本实施例的DEM数据来源于NASA的STRM数据，但数据来源不局限于STRM，也可以采用ASTER GDEM或其他。

1.3.获取覆盖巴林右旗区域的气象数据，包括气温、降水、风速、风向等气象要素数据；上述气象数据分别采用克里金插值方法进行插值，生成各气象要素栅格数据，并利用巴林右旗矢量边界对气象栅格数据进行裁剪。

1.4.获取覆盖巴林右旗区域的土地覆盖类型数据(实施例的数据来源于2014年)，首先进行配准，再利用巴林右旗矢量边界进行裁剪，获取巴林右旗的土地覆盖类型图，该土地覆盖类型图包括植被类型、建设用地类型等属性。

1.5.获取覆盖巴林右旗区域的遥感图像，选取红外波段与近红外波段计算植被归一化指数(NDVI)，并绘制NDVI图；在此基础上，运用像元分解模型，获取巴林右旗的植被覆盖度数据，经过裁剪，得到巴林右旗的植被覆盖度图。

其中，NDVI的计算公式为式中NIR，RED为近红外和红光波段反射率。

像元分解模型为其中，NDVI_soil为裸土或无植被覆盖区域的NDVI值，即无植被像元的NDVI值；而NDVI_veg则代表完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值。

需要说明的是，根据巴林右旗区域的气候和植被特征，选取植物生长旺季7-8月的遥感影像，再根据影像清晰程度进行筛选，最终选取了2014年7-8月美国Landsat卫星TM7遥感影像。但是，本发明并不局限于Landsat数据，也不局限于上述时相，具体选取的遥感影取决于研究区的气候和植被特征，以及遥感数据的可获取性以及遥感影像的质量。

2.计算生态功能要素，包括水源涵养、防风固沙、土壤保持和生物多样性维持这四个方面，分别通过各自的公式进行计算，并分别生成栅格图。计算步骤利用栅格计算器进行，具体利用ArcGIS软件中的“raster calculator”功能，其中计算公式中的植被盖度、坡度、降雨量等参数数据均来自步骤1。

2.1.水源涵养量的计算。

水源涵养量的计算采用降水贮存量法，主要是以森林和草地生态系统的蓄水效应来表征水源涵养功能，计算公式为Q＝A·J·R，最终生成表征水源涵养量的图。

式中，J＝J₀·K，R＝R₀-R_g，R_g＝-0.3187×P_g+0.36403；

Q为与裸地相比较，不同生态系统涵养水分的增加量(mm/(hm²·a^-1))；

A为生态系统面积(hm²)，基于步骤1.1中的矢量边界计算；

J为研究区多年均产流降雨量(P>20mm)(mm)；

J₀为研究区多年均降雨总量(mm)，基于步骤1.3中的降水栅格数据计算；

K为研究区产流降雨量占降雨总量的比例，该数值可以预先通过科研成果或者经验获得。本实施例采用赵同谦(2004)等的研究结果，以秦岭-淮河一线为界限将全国划分为北方区和南方区，北方区K取0.4，南方区K取0.6；

P_g为植被覆盖度，采用步骤1.5的计算结果；

R₀为产流降雨条件下裸地降雨径流率，即当植被覆盖度P_g＝0时的R_g值，依据公式R_g＝-0.3187×P_g+0.36403，可得出R₀值为0.36403；

R_g为产流降雨条件下生态系统降雨径流率，依公式R_g＝-0.3187×P_g+0.36403获得；

R为与裸地(或皆伐迹地)比较，生态系统减少径流的效益系数，依据公式R＝R₀-R_g计算得到。

2.2.防风固沙量的计算。

防风固沙量的计算由潜在土壤侵蚀量和现实土壤侵蚀量估算。国标土壤侵蚀分类分级标准SL190-2007规定了风蚀区土壤侵蚀强度的分级标准，并根据植被覆盖度给出了土壤侵蚀强度的分级指标。巴林右旗研究区的风蚀土壤侵蚀量主要受植被覆盖度影响，根据表1确定不同植被覆盖度下的防风固沙量，基于植被覆盖度图，绘制巴林右旗防风固沙量的空间分布图，生成表征防风固沙量的图。

表1防风固沙量确定单位：t··km^-2·a^-1

2.3土壤保持生态功能的计算。

土壤保持功能的计算是根据当地气候特点(降雨)、地形特点(坡度坡长)、植被特点(覆盖因子)、土壤类型(土壤可蚀性因子)和土地管理水平等因素对生态系统土壤保持功能的强弱进行评估。

土壤保持功能一般采用通用水土流失方程(USLE)进行评价。在计算时，需要根据不同采样点实测的数据对模型模拟结果进行验证，并且进一步修正参数。最终生成表征土壤保持功能的图。

USLE_x＝R_x·K_x·LS_x·C_x·P_x

式中：USLE_x表示栅格x的土壤侵蚀量，R_x为降雨侵蚀力，K_x为土壤可蚀性，LS_x为坡度-坡长因子(已通过步骤1.2获得)，C_x为植被覆盖度(通过步骤1.5获得)，P_x为管理因子。

其中，降雨侵蚀力R_x由Fouriner指数计算，计算公式为上式中，R降雨侵蚀力、j表示月降水、J表示年降水、i表示月份。

土壤可蚀性因子K_x采用土壤风蚀方程中土壤可蚀性因子计算模型$K_x=\frac{29.09+0.31sa+0.17si+0.33(sa/cl)-2.59OM-0.95{\mathrm{CaCO}}_3}{100}$

式中，K_x表示土壤可蚀因子；sa表示土壤粗砂含量；si表示土壤粉砂含量；cl表示土壤粘粒含量；OM表示有机质含量；CaCO₃表示碳酸钙含量。土壤颗粒含量(sa、si、cl)和有机质含量数据来源于中国土壤特征数据集，由中国西部环境与生态科学数据中心(http://westdc.westgis.ac.cn)提供；土壤碳酸钙含量来自于全国1:400万土壤碳酸钙含量分布图，由地球系统科学数据共享服务网(http://www.geodata.cn)提供。管理因子P根据土壤利用类型，通过表2进行确定。

表2局部区域P值

2.4生物多样性维持功能的计算

生物多样多样性维持功能可以采用Invest-biodiversity模型进行计算，生物多样性维持功能用生境质量好坏代表生物多样性的持续性、恢复能力，因此模型中采用生境质量指数进行评价，最终生成表征生物多样多样性维持功能的图。

Invest-biodiversity评价模型中输入土地覆盖类型图、威胁因子、威胁因子图层、地类对于威胁因子的敏感度图层计算得到生境质量指数结果，评价模块运算以栅格数据作为评价单元，所有栅格大小为30m分辨率。

其中，土地覆盖类型的分类参照全国生态环境十年变化遥感调查与评估项目，土地类型分类，包括6大类共21种小类土地利用类型，包括森林(落叶阔叶林、常绿针叶林、落叶阔叶灌木林等)，草地(草甸、草丛、草本、绿地等)，湿地(河流、湖泊、水库、坑塘、运河水渠等)，农田(旱地、水田等)，建设用地(居住地、厂矿、交通用地等)，未利用地(裸岩、裸土、稀疏林和稀疏草地等)。

其中，威胁因子为对生物多样性造成威胁的图层，包括耕地、城镇、工业用地、铁路、高速公路、国道和省道等因子，其中除土地覆盖类型因子外，道路数据来自国家测绘地理信息局，威胁因子参照invest2.6.5Beta User’s Guide中研究成果定义危险因子的最大影响距离、权重及衰退现行相关指数如图1，在此基础上形成威胁因子图。

其中，地类对于威胁因子的敏感度图层主要是依据生态学和景观生态学的基本理论及保护生物多样性的基本原则来确定，参照invest2.6.5Beta User’s Guide中生态威胁因子划分标准，在生物多样性模型中为各土地利用类型对于威胁因子敏感度赋值，取值范围为0-1，如森林、草地和湿地等天然土地利用类型敏感度较高，耕地、建设用地等人工土地利用类型敏感度较低，裸土裸岩等未利用地敏感度最低，对于各土地利用类型自然属性一列的赋值是将天然地类和人工地类分别赋值为1、0，见图2，在此基础上形成地类对于威胁因子的敏感度图层。

3.判定生态退化程度，生成研究区的生态退化程度图。

依据理论NDVI与实际NDVI比值，参照退化等级划分标准(表3)，确定研究区的生态退化等级，生成生态退化等级图。

表3退化等级划分标准

其中，理论NDVI的确定方法如下：获取1985-1990年间，每年最大降水量对应的月份；对研究区域进行随机取样，实施例中共取样点1400个；然后在1985-1990年间，每年出现最大降水量的月份的NDVI图像中，将上述采样点对应位置的NDVI值提取出来，再通过对NDVI值作平均值处理，求取代表研究区域、最大降水量月份的平均NDVI值；将最大降水量与平均NDVI值分别作为X、Y轴作出散点图，分析降水量和NDVI的散点关系，采用统计分析软件(如SPSS)对其进行线性回归分析，建立回归模型；将2014年的最大月降水量作为自变量分别代入回归模型，求得该降水量下的理论NDVI植被指数，该理论NDVI表示在该降水量(气候条件)下，植被生长所能达到的最佳状态。其中，选择1985-1990年作为求取理论NDVI的基础，是因为该时间段中研究区的人为干扰最小。但是，本发明并不局限于这个时间段，其他实施例应根据研究区的实际情况和遥感数据的可获取状况进行选择。NDVI的计算是基于步骤1.5中公开的方法，计算NDVI的遥感数据源来源于TM7。

实际NDVI依据遥感影像进行计算，实际NDVI已在步骤1.5中得到结果。

4.划定各级生态修复空间分区(如图3)，并分别进行赋值。

生态修复分区主要受到3个因素的影响，包括立地条件(气候气象、地形地貌、植被覆被)、主导生态功能变化、生态退化等级。

4.1确定研究区生态修复空间分区的影响要素(图4)。

基于不同空间尺度、不同气候条件、不同下垫面，选取不同的空间分区指标要素。实施例的巴林右旗面积约1万平方公里，气候上属于干旱、半干旱气候区，气候气象因素在空间分布上的差异不明显，因此未考虑气候因素。实施例中确定的生态修复空间分区的影响要素，包括立地条件(仅考虑地形地貌、植被覆被)、主导生态功能变化、生态退化等级。需要说明的是，本发明的其他实施例根据具体的地理位置和气候条件，可以选取具体应用的影响要素。

4.2划分一级区。

一级区为高级区划单位，区域内存在共同的环境特征和生态问题，土地资源主导利用方向基本一致，而且区域特征、分布范围比较明显，具有空间尺度的完整性和连续性。就巴林右旗为例，选用立地条件中的海拔和土地覆盖类型数据中森林、草地和灌丛的属性值作为一级区的分级指标，并生成生态修复一级区图。巴林右旗区域的海拔变化范围从355m-1900m，以1000米-1900米范围内地貌为中山区，低于1000米为丘陵区，低于500m为平原区。根据上述规则进行划分，巴林右旗地区的北部为山地区，中部为丘陵区，南部为平原区。中山区主要植被类型为森林、丘陵区主要植被类型为灌丛，平原区主要植被类型为草地。

实施例将DEM作为地形地貌中海拔的数据源，应用空间分析技术，依据上述界定的海拔范围进行划定并赋值，如此初步划定了一级区，并生成生态修复一级区图。具体可以利用ArcGIS中的spatial analyst\reclassify功能。本实施例中一级区的赋值为4位数，例如可以是1000、2000、3000等。

4.2划分二级区。

二级区为中级区划单位，也是一级区单位的细分，主要反映中尺度下生态修复空间生态功能的分异规律，因此选取水源涵养、防风固沙、土壤保持和生物多样性这4类生态功能要素中的主导生态功能要素作为划分指标，生成主导生态功能空间分布图。

主导生态功能要素通过以下方法确定：

步骤2中通过各自的公式，对巴林右旗水源涵养、防风固沙、土壤保持和生物多样性这四大生态功能要素分别进行了计算，并分别生成了栅格图A、B、C、D。由于这4个生态功能单位不一致，无法进行比较和栅格运算，因此采用极差标准化法分别将4个栅格数据进行标准化，使得4个栅格数据的取值在0-1的标准化值范围内。标准化公式如下：

标准化值

X_min、X_max分别表示指标值中的最小值和最大值)；指标其他数值均在最大值和最小值之间。

将经过标准化的图层A、B、C、D分别进行栅格属性值重分类，分为以下4个等级，即一般重要(0-0.3)、中等重要(0.3-0.6)、重要(0.6-0.8)、极重要(0.8-1)，生成A1,B1,C1,D1这4个新栅格图；再将A1、B1、C1、D1这4个栅格数据进行叠加分析，将各个栅格对应位置上，A1～D1中最大等级对应的生态功能作为主导生态功能要素，并通过栅格属性值重分类进行对应赋值，将具有相同赋值的栅格合并为一个斑块，生成新的主导生态功能空间分布图。其中，本实施例中二级区的赋值为3位数，例如防风固沙功能“100”、土壤保持功能“200”、水源涵养功能“300”、生物多样性维持功能“400”。

上述栅格属性值重分类可以采用ArcGIS软件中的空间分析模块中模块中的reclassfy功能；叠加分析也可以在ArcGIS软件中实现。

4.3划分三级区。

按照步骤3中对生态退化程度的判定，已经将巴林右旗地区划定为未退化、轻度退化、中度退化、重度退化这4种退化等级。其中，在未退化区进一步采用自然断裂法，按照实际NDVI值大小分为3类，依据NDVI数值从大到小将未退化区进一步划分为优先保护区、重点保护区和一般保护区。此外，将步骤3中确定的轻度退化划定为预防修复区，中度退化划定为调控修复区，重度退化划定为重点修复区，将上述新划定的区域进行赋值，将具有相同赋值的栅格合并为一个斑块，最终生成生态修复三级区图。

其中，本实施例中三级区的赋值为2位数，例如各分区的赋值可以是优先保护“10”、重点保护“20”、一般保护区“30”、预防修复区“40”、调控修复区“50”、重点修复区“60”，通过栅格属性值重分类生成生态修复三级区图。上述自然断裂法可以采用ArcGIS软件中的natural breaks功能。

5.生成生态修复分区图并命名。

将步骤4中生成的生态修复一级分区图、二级分区图、三级分区图进行叠加运算，生成生态修复分区图。上述叠加运算是针对各分区图中赋值的code层。实施例最终生成的生态修复分区图的code统一为四位数，末尾数为相同立地条件、相同主导功能、相同退化程度下栅格的顺序数，生态修复分区命名依据该三级划分的四位数命名，如：①北部高山森林水源涵养生态修复预防小区(1341)；

②中部丘陵草地防风固沙生态修复重点小区(2161)；

③南部平原开发用地水土保持生态修复调控小区(3252)。

上述叠加运算可以采用ArcGIS软件中的raster caculator命令。巴林右旗生态修复分区结果见表4，空间分布图见附图5。

表4巴林右旗生态功能保护与恢复区划表

综上所述，本发明实施例提供了一种结合主导生态功能与生态退化程度的生态修复分区方法，本发明在遵从研究区立地条件的前提下，以提升生态功能为目标，并结合生态退化程度判定生态修复难易程度，构建生态修复空间分区指标要素，进而划定生态修复三级分区，克服了以往生态功能区划只重视生态系统主要服务功能空间布局的局限性，提高了未来生态修复工程的针对性和精确性，节省了人力、物力。

本领域技术人员可以理解，实现上述实施例方法的全部或部分流程，可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于计算机可读存储介质中。其中，所述计算机可读存储介质为磁盘、光盘、只读存储记忆体或随机存储记忆体等。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。