CN116757357B - 一种耦合多源遥感信息的土地生态状况评估方法 - Google Patents

Abstract

本发明公开了一种耦合多源遥感信息的土地生态状况评估方法，包括从“生境状况‑结构状况‑胁迫状况‑土地利用效益”角度构建多层级综合评价指标体系；确定评价单元，充分利用多源遥感数据提取并优化评价指标；对指标数据进行标准化处理；采用主客观综合权重法测算评价指标的权重；利用综合指数法评估土地生态状况，进行等级划分和分析。本发明利用多源遥感数据，更客观地表征了土地生态状况的自然属性；指标的优化提高了土地生态状况评估的客观性、科学性和合理性；综合考虑土地的生物环境、结构和弹性、外部胁迫及提供效益的能力等因素，构建了一套科学客观的土地生态状况评价指标体系，过程明晰，具有较强的普适性，易于在区域尺度中推广应用。

Description

一种耦合多源遥感信息的土地生态状况评估方法

技术领域

本发明涉及一种耦合多源遥感信息的土地生态状况评估方法，属于生态建设、生态环境保护、生态承载能力评价、国土资源安全、国土空间规划等应用领域。

背景技术

土地资源是人类赖以生存的物质基础和社会发展的核心，其生态状况对人类生存和发展至关重要。随着社会经济的飞速发展，人类对土地资源的利用日益加剧，土地资源的高强度利用为社会发展带来了巨大的经济效益，同时水土流失、土地污染与退化、土地沙漠化等土地生态环境问题也越来越严重。因此，科学客观评价土地生态状况，对构建高质量发展的国土空间布局，实现土地资源可持续发展具有重要的现实意义。

目前，土地生态状况研究主要有土地生态安全、生态质量、生态敏感性、生态脆弱性、土地退化和土地生态预测等众多内容，由于土地生态状况研究对象的广泛性、内容的丰富性和尺度的多样性，至今尚未形成一条较为系统的土地生态状况评价体系。目前最常用的指标体系模型为生态基础-生态结构-生态胁迫-生态效益模式，但该模式只考虑了土地的生态指标，忽略了土地利用的社会经济活动。

其次，已有研究评价指标体系中的指标有较大差异，涉及站点数据、统计数据及遥感影像等多种数据。大多研究中的水热条件主要使用气象站点数据，但站点密度、分布状况以及下垫面复杂程度会造成空间插值数据具有很大的不确定性；统计数据是以行政区划为统计单元，在应用范围上具有一定程度的局限性，空间特性较隐含，而且无法反应区域内部的空间差异和分布特征。而且某些评价指标不尽合理，例如温度是反映水热状况的关键参数，大多数研究使用站点年均温来表示，未协同考虑温度强度和温度作用时间对植物生长、发育的影响；土地的利用结构大多使用单独的林地、草地等多种用地类型的占比来表示，能较好地反映土地利用中某一侧面的特点，但不利于反映土地利用的总体结构，而且不同区域之间也无法进行横向比较。

从指标赋权和评价方法模型上看，目前处于不断探索阶段。指标赋权主要包含主观、客观、综合三种，主观赋权法的主观性强，容易忽略掉某些信息；客观赋权法则基于数据反映的信息水平进行赋权，会产生违背社会现实的权重分配，而综合赋权法则综合利用主观、客观的方法，吸收了两者优点，同时弥补了缺点。较为常用的评价方法模型包括数学模型法、空间模型法和景观生态模型法，其中数学模型法中的综合指数法是应用最多的方法，该方法的优点是过程简单、不损失指标信息，评价指标横向与纵向的对比分析较简便。

发明内容

针对现有技术方法存在的不足，本发明提出了一种耦合多源遥感信息的土地生态状况评估方法，该方法综合考虑了土地生态的生物环境、结构和弹性、外部胁迫及提供效益的能力等因素，构建了一套科学客观的土地生态状况评价指标体系；并充分利用多源遥感数据提取和优化评价指标，提高了土地生态状况评估的客观性、科学性和合理性；最后利用主客观综合权重法和综合指标法科学全面的评估了土地生态状况，过程明晰、具有较强的普适性，易于在区域尺度中推广应用。

为解决上述问题，本发明所采取的技术方案是：

一种耦合多源遥感信息的土地生态状况评估方法，包括以下步骤：

S1，基于17个评价指标构建多层级综合评价指标体系；

S2，采用1km×1km格网对研究区进行评价单元的划分；

S3，利用多源遥感数据并采用步骤S2中的1km×1km格网对步骤S1中的17个评价指标进行提取和优化，得到17个评价指标初始数据；

S4，对步骤S3中的评价指标初始数据进行极差标准化处理，得到17个标准化后的评价指标极标数据，极差标准化的计算公式如下：

若评价指标为正向：

若评价指标为逆向：

式中，x_maxi、x_mini分别为第i个评价指标初始数据的最大值和最小值，x_mi为第i个评价指标初始数据在第m个评价单元的值，y_mi为第i个评价指标极标数据在第m个评价单元的值；

S5，针对步骤S4中标准化后的评价指标极标数据，利用层次分析法和熵权法主客观综合赋权测算各评价指标极标数据的权重；

S6，针对步骤S4中标准化后的评价指标极标数据和步骤S5中各评价指标极标数据的权重，利用综合指数法评估土地生态状况，得到土地生态状况综合指数LESI，LESI计算公式如下：

式中，LESI为土地生态状况综合指数；X_ij为第i个评价单元的第j个评价指标极标数据的值；w_j为评价指标极标数据X_ij的权重；z为评价指标极标数据的数量；

S7，针对步骤S6中的土地生态状况综合指数LESI，利用自然段点法划分土地生态状况等级，得到土地生态状况等级空间分布结果。

优选的，步骤S1中多层级综合评价指标体系包括生境状况、结构状况、胁迫状况和土地利用效益4个准则层，4个准则层包括10个要素层，10个要素层包括17个评价指标。

优选的，所述生境状况包括立地条件、土壤条件和水热条件；所述结构状况包括土地利用结构和生态景观；所述胁迫状况包括人口压力和人为干扰；所述土地利用效益包括社会效益、经济效益和生态效益；

其中，所述立地条件为地形位指数，所述土壤条件为土壤质地，所述水热条件包括年降水量和活动积温，所述土地利用结构为土地利用综合指数，所述生态景观包括景观多样性指数、景观蔓延度指数和河网密度，所述人口压力为人口密度，所述人为干扰包括路网密度和不透水盖度，所述社会效益为地均粮食产量，所述经济效益为地均GDP，所述生态效益包括人均生态用地、生态服务价值、植被覆盖度和植被总初级生产力。

优选的，步骤S1中的17个评价指标分别为：地形位指数、土壤质地、年降水量、活动积温、土地利用综合指数、景观多样性指数、景观蔓延度指数、河网密度、人口密度、路网密度、不透水盖度、地均粮食产量、地均GDP、人均生态用地、生态服务价值、植被覆盖度、植被总初级生产力。

优选的，步骤S2中1km×1km格网评价单元是先以MODIS LST数据为环境变量，利用ArcGIS Create Fishnet工具生成格网；再利用研究区矢量进行裁剪得到。

优选的，步骤S3中17个评价指标的提取和优化如下：

(1)地形位指数，利用规范化处理后的DEM遥感数据和坡度计算得到研究区内的地形位指数初始数据，计算公式如下：

式中，T为地形位指数；E和S分别表示区域范围内任意一点M的高程值和坡度值；和分别表示M点所在区域的平均高程值和平均坡度值，其中，坡度值由DEM数据通过ArcGISSlope工具计算得到；

(2)土壤质地，使用土壤粉砂含量，经统一的规范化处理后得到研究区内的土壤质地初始数据；

(3)年降水量，基于卫星遥感降水数据，对同年度区域所有日降水数据进行求和，经统一的规范化处理后得到研究区内的年降水量初始数据；

(4)活动积温，通过NDVI时间序列数据确定的植被生长开始日期和结束日期，计算植被生长开始日期至结束日期的地表温度LST累计和，经统一的规范化处理后得到研究区内的活动积温初始数据；

(5)土地利用综合指数，基于遥感解译的土地利用数据，从三生空间的角度出发，构建一种土地利用综合指数模型，并通过一种土地利用综合指数模型计算得到土地利用综合指数初始数据；

(6)景观多样性，使用Shannon Weaver指数度量，提取步骤为：利用最优尺度景观单元格网对遥感解译的土地利用数据进行测算，经统一的规范化处理得到景观多样性初始数据，计算公式如下：

式中，SHDI为景观多样性指数，p_i为景观斑块类型i所占的面积百分比，q为景观中的斑块类型总数目；

(7)景观蔓延度，提取步骤为：利用最优尺度景观单元格网对遥感解译的土地利用数据进行测算，经统一的规范化处理得到景观蔓延度初始数据，计算公式如下：

式中，CONTAG为景观蔓延度指数，p_i为景观斑块类型i所占的面积百分比，g_ik为i类型斑块和k类型斑块毗邻的数目，m为景观中的斑块类型总数目；

(8)河网密度，指1km×1km评价单元内河流的长度，提取步骤为：首先提取遥感解译的面状河流数据的中心线，然后利用1km×1km的网格评价单元统计每个网格内的河流长度，再转化成栅格数据，经统一的规范化处理后得到河网密度初始数据；

(9)人口密度，在分县人口统计数据的基础上，综合考虑了与人口密切相关的土地利用类型、夜间灯光亮度、居民点密度因素，将以行政区为基本统计单元的人口数据展布到空间格网上，从而实现人口的空间化，经统一的规范化处理后得到人口密度初始数据；

(10)路网密度，指1km×1km评价单元内所有道路的总长度；提取步骤为：首先提取遥感解译的县道、省道、国道和铁路等面状道路数据的中心线，再利用1km×1km的网格统计每个网格内的道路总长度，最后转化成栅格数据，经统一的规范化处理后得到路网密度初始数据；

(11)不透水盖度，指1km×1km评价单元内不透水面的面积，提取步骤为：首先基于遥感解译的不透水数据，利用1km×1km的网格统计每个网格内不透水面面积，再转化成栅格数据，经统一的规范化处理后得到不透水盖度初始数据；

(12)地均粮食产量，在粮食产量分县统计数据的基础上，构建一种基于MODIS GPP数据实现粮食产量空间化的模型，并通过一种基于MODIS GPP数据实现粮食产量空间化的模型计算得到地均粮食产量初始数据；

(13)地均GDP，在分县GDP统计数据的基础上，综合考虑与人类经济活动密切相关的土地利用类型、夜间灯光亮度、居民点密度因素，将以行政区为基本统计单元的GDP数据展布到栅格单元上，从而实现GDP的空间化，经统一的规范化处理后得到地均GDP初始数据；

(14)人均生态用地面积，指1km×1km评价单元内生态用地面积与人口的比值，提取步骤为：首先基于遥感解译的生态用地，利用1km×1km网格统计每个格网内土地生态用面积，再转化成栅格数据，最后经统一的规范化处理后再除以人口密度数据，得到人均生态用地初始数据；

(15)生态服务价值，提取步骤为：首先计算1km×1km评价单元内不同地类面积和生态服务价值系数乘积的和，再经统一的规范化处理后得到生态服务价值初始数据，计算公式如下：

式中，ESV为生态系统服务价值总量，单位：万元；A_s是研究区第s种地类分布面积；VC_s为第s种地类的生态服务价值系数；

(16)植被覆盖度，提取步骤为：首先对同年度区域所有遥感影像进行精确的辐射定标、大气校正和几何配准；然后利用云层指数CCI和阴影指数YYI消除云和阴影的影响，计算公式如下：

CCI＝(ρ_R+ρ_G+ρ_B)×ρ_NIR；

式中，ρ_R、ρ_G、ρ_B、ρ_NIR、ρ_SWIR分别为红波段、绿波段、蓝波段、近红外波段和短波红外波段的光谱反射率值，NDVI为归一化植被指数；

最后利用最大值合成法提取年度内所有NDVI数据的最大值，并采用像元二分法计算植被覆盖度，经统一的规范化处理后得到植被覆盖度初始数据，计算公式如下：

式中，VFC为像元的植被覆盖度，NDVI为该像元上的归一化植被指数，NDVI_soil为无植被覆盖像元的NDVI值，NDVI_veg为完全被植被覆盖的像元NDVI值，其中，NDVI_veg和NDVI_soil分别为NDVI累积百分数5％和95％对应的NDVI值；

(17)植被总初级生产力，提取步骤为：首先对MODIS GPP时间序列数据进行线性插值和S-G滤波；然后计算年均值，经统一的规范化处理后得到研究区内的植被总初级生产力初始数据；

其中统一的规范化处理为转投影、重采样和裁剪，以使17个评价指标初始数据的空间信息与1km×1km格网评价单元一致。

优选的，一种土地利用综合指数模型的构建包括以下步骤：

首先，从三生空间角度将遥感解译的土地利用数据分类为生态空间、生产空间和生活空间，其中生产空间是以提供工业品、农产品和服务产品为主导功能的区域；生活空间是以提供人类居住、消费、休闲和娱乐为主导功能的区域；生态空间是以提供生态产品和生态服务为主导功能的区域；

其次，对三生空间进行分级赋值，生态空间赋值3，生产空间赋值2，生活空间赋值1；

最后，利用1km×1km网格统计每个格网内各用地空间类型等级值和面积百分比乘积的和，并转化成栅格数据，得到土地利用综合指数初始数据，计算公式如下：

式中，L_c为土地利用综合指数，A_i为第i种用地空间类型的分级赋值，C_i为第i种用地空间类型的面积百分比，n为用地空间类型的总数量。

优选的，景观多样性和经蔓延度中的最优尺度景观单元格网确定方法如下：

首先，利用ArcGIS Create Fishnet创建1km、2km、2+N…50km共26个窗口尺度的景观单元网格；

其中，N为大于等于2且小于等于46的偶数；

其次，分别利用26个窗口尺度的景观单元网格裁剪土地利用数据，得到相应的景观单元的土地利用类型数据；

然后，将景观单元土地利用数据批量输入Fragstats4.2软件计算景观多样性指数SHDI；

最后，通过尺度效应曲线选取最优尺度的景观单元格网，用于提取景观多样性和景观经蔓延度。

优选的，一种基于MODIS GPP数据实现粮食产量空间化的模型的构建包括以下步骤：

首先，对全年GPP时间序列数据进行线性插值和S-G滤波，以及统一的规范化处理；

其次，提取每个耕地像元的全年GPP总和，记为GPP_T；

然后，利用县域行政区矢量统计各县范围内所有耕地像元的GPP_T的总和；

最后，利用像元GPP_T占县域GPP_T的总和的比，再乘以县域粮食产量，实现粮食产量的空间化，得到地均粮食产量初始数据，计算公式如下：

式中，Y_ik为第k个县域行政区第i个像元的粮食产量，GPP_T_ik为第k个县域行政区第i个像元的GPP年总和值，u为第k个县域行政区的像元总数，Y_k为第k个县域行政区粮食产量。

本发明的有益效果在于：

1)本发明的优点是综合考虑了土地的生物环境、结构和弹性、外部胁迫状况及提供效益的能力等因素，构建了一套“生境状况-结构状况-胁迫状况-土地利用效益”多层级综合评价指标体系；

2)建立了一套包含17个因子的评价指标体系，此外还对活动积温、土地利用结构、景观多样性、景观蔓延度和地均粮食产量等指标的计算方法进行了优化，提高了土地生态状况评估的客观性、科学性和合理性。

3)相比现有技术方法，本发明充分利用多源遥感数据，改变了以往过多依赖于站点数据而对土地生态状况自然属性重视不够的不足之处；数据源除土壤质地及人口粮食产量和GDP统计数据外，均为遥感数据能更客观的反应生态状况，且易获取。

4)提出了一种结合遥感数据的社会经济统计数据空间化处理方法，能更好的反应区域内部的空间差异和分布特征，并有效的解决了统计数据与遥感数据空间不匹配的问题。

附图说明

图1为本申请实施例的流程示意图；

图2为本申请实施例中构建的土地生态状况多层级综合评价指标体系；

图3为本申请实施例中1km×1km格网评价单元；

图4为本申请实施例中景观多样性和景观蔓延度尺度效应曲线。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

本实施例以河北省某年度耦合多源遥感数据的土地生态状况评估为例进行说明。多源遥感数据包括DEM、GSMAP卫星遥感、MODIS遥感、Landsat遥感，这些均为官方网站数据；植被覆盖度遥感反演数据，土地利用遥感解译数据、路网遥感解译数据和不透水遥感解译数据的精度均大于85％，可以满足实际应用要求。

参照图1，本发明公开了一种耦合多源遥感信息的土地生态状况评估方法，其包括以下步骤：

S1，在充分认识土地生态评价特点和内涵的基础上，从“生境状况-结构状况-胁迫状况-土地利用效益”角度，并基于17个评价指标构建多层级综合评价指标体系；

具体的，从土地生态系统本身出发，土地生态状况为是土地的生物环境、结构和弹性、外部胁迫状况以及提供效益的能力等因素综合作用的结果。一个良好的土地生态状况，应当是具有较好的生物环境，较稳定的结构状况，受到的外部胁迫状况较少，同时能带来较大的效益。

多层级综合评价指标体系包括4个准则层，4个准则层包括10个要素层，10个要素层包括17个评价指标。

其中，4个准则层分别为生境状况、结构状况、胁迫状况和土地利用效益；

其中，生境状况包括立地条件、土壤条件和水热条件；结构状况包括土地利用结构和生态景观；胁迫状况包括人口压力和人为干扰；土地利用效益包括社会效益、经济效益和生态效益；10个要素层分别为立地条件、土壤条件、水热条件、土地利用结构、生态景观、人口压力、人为干扰、社会效益、经济效益和生态效益。

其中，立地条件为地形位指数，土壤条件为土壤质地，水热条件包括年降水量和活动积温，土地利用结构为土地利用综合指数，生态景观包括景观多样性指数、景观蔓延度指数和河网密度，人口压力为人口密度，人为干扰包括路网密度和不透水盖度，社会效益为地均粮食产量，经济效益为地均GDP，生态效益包括人均生态用地、生态服务价值、植被覆盖度和植被总初级生产力。17个评价指标为：地形位指数；土壤质地；年降水量；活动积温；土地利用综合指数；景观多样性指数；景观蔓延度指数；河网密度；人口密度；路网密度；不透水盖度；地均粮食产量；地均GDP；人均生态用地；生态服务价值；植被覆盖度；植被总初级生产力。

其中，生境状况是土地生态系统功能正常运行的基础，结构状况反应了土地对外界干扰的抵抗能力，胁迫状况反应了外部环境对土地生态的压力和胁迫，土地利用效益则是土地生态可持续发展的动力和支撑。详细见图2和表1。

表1

其中，“+”代表正向指标，“-”代表逆向指标。

S2，为了保证单元生态信息完整的同时提升评价精度，采用1km×1km格网对研究区进行评价单元的划分；

具体的，在空间尺度为市级或更大尺度时，常以行政区作为尺度单元，但这种方法过于笼统和粗糙，淹没了小尺度的大量细节信息，以1km×1km格网划分评价单元可以保证单元生态信息完整的同时提升评价精度。

操作步骤为：首先以1km空间分辨率MODIS LST数据为环境变量，利用ArcGISCreate Fishnet工具生成1km×1km格网；再利用河北省矢量范围进行裁剪，得到1km×1km格网评价单元，见图3。

S3，充分利用多源遥感数据，采用步骤S2中的1km×1km格网对步骤S1中的17个评价指标进行提取和优化，得到17个评价指标初始数据；

具体的，各评价指标初始数据来源见表1，提取过程中均需对17个评价指标进行统一的规范化处理，统一的规范化处理包括转投影、重采样和裁剪等工作，保证17个评价指标初始数据的空间信息和范围与1km×1km格网评价单元一致。计算方法如下：

(1)地形位指数，利用ArcGIS栅格计算器对规范化处理后的DEM遥感数据和坡度进行计算，得到研究区内的地形位指数初始数据。计算公式如下：

式中，T为地形位指数；E和S分别表示区域范围内任意一点的高程和坡度值；和分别表示该点所在区域的平均高程和平均坡度值。其中，坡度由DEM数据通过ArcGIS Slope工具计算得到。

具体的，地形差异是至关重要的自然本底要素，直接影响区域的气候条件和人类对土地的利用方式；而高程和坡度，一定程度上决定着人类的农业生产和城市的建设活动，是区域地形条件最主要的表征。

(2)土壤质地，使用土壤粉砂含量。经统一的规范化处理后得到研究区内的土壤质地初始数据。

具体的，土壤粉砂含量与土壤通气、保肥、保水状况及耕作的难易有密切关系，是土地生态状况最基础的指标。

(3)年降水量，基于GSMAP卫星遥感降水数据，对2020年区域所有日降水数据进行求和，经统一的规范化处理后得到研究区内的年降水量初始数据。

具体的，人类活动及自然界的许多变化，同土地资源的作用最直接最集中的体现在水、热两个方面。地面气象站点观测是获取降水数据的常规方法，但受限于站点密度、分布状况以及下垫面复杂程度的影响，通过空间插值生成空间连续分布的降水数据具有很大的不确定性；卫星遥感降水数据以其高时空分辨率优势特征，逐渐成为重要数据源。

(4)活动积温，主要通过NDVI时间序列数据确定的植被生长开始日期和结束日期，计算植被生长开始日期至结束日期的地表温度LST累计和，经统一的规范化处理后得到研究区内的活动积温初始数据。

具体的，温度是反映水热条件的另一个关键参数。目前大多数研究使用地面气象站点年均温来表示，空间插值数据的不确定性会对结果产生一定的影响，因此使用遥感数据反演的地表温度来代表区域的热量条件。

(5)土地利用综合指数，基于遥感解译的土地利用数据，从三生空间的角度出发，构建一种土地利用综合指数模型，通过土地利用综合指数模型计算得到土地利用综合指数初始数据。

具体的，目前常利用多个指数来反映土地利用的方式和结构，能较好地反映土地利用中某一侧面的特点，但不利于反映土地利用的总体结构，而且不同区域之间也无法进行横向比较从而进行其区域分异的分析。因此构建土地利用综合指数定量的表达土地利用结构，构建方法如下：

首先，从三生空间角度将遥感解译的土地利用数据分类为生态空间、生产空间和生活空间，其中生态空间包括林地、草地、水域、未利用土地；生产空间为耕地；生活空间为城乡、工矿、居民用地；

其次，对三生空间进行分级赋值，生态空间赋值3，生产空间赋值2，生活空间赋值1；

最后，利用1km×1km网格统计每个格网内各用地空间类型等级值和面积百分比乘积的和，并转化成栅格数据，得到土地利用综合指数初始数据。

计算公式如下：

其中，L_c为土地利用综合指数，A_i为第i种用地空间类型的分级赋值，C_i为第i种用地空间类型的面积百分比，n为用地空间类型的总数量。

(6)景观多样性，使用Shannon Weaver指数度量。提取步骤为利用最优尺度景观单元格网对遥感解译的土地利用数据进行测算，经统一的规范化处理得到景观多样性初始数据。计算公式如下：

式中，SHDI为景观多样性指数，p_i为景观斑块类型i所占的面积百分比，q为景观中的斑块类型总数目。

具体的，景观要素具有尺度效应，小尺度条件无法体现大的景观生态背景，大尺度背景又可能会忽略对细节的描述，因此需要确定最优尺度，详细的确定方法如下：

首先利用ArcGIS Create Fishnet创建1km、2km、2+N…50km共26个窗口尺度的景观单元网格，其中，N为大于等于2且小于等于46的偶数；其次分别利用26个窗口尺度的景观单元网格裁剪土地利用数据，得到相应的景观单元的土地利用类型数据；然后将景观单元土地利用数据批量输入Fragstats4.2软件计算景观多样性指数SHDI；最后通过尺度效应曲线选取最优尺度的景观单元格网为10km，见图4。

(7)景观蔓延度，提取步骤为利用最优尺度景观单元格网对遥感解译的土地利用数据进行测算，其中，最优尺度景观单元格网可选用10km，经统一的规范化处理得到景观蔓延度初始数据。计算公式如下：

式中，CONTAG为景观蔓延度指数，p_i为景观斑块类型i所占的面积百分比，g_ik为i类型斑块和k类型斑块毗邻的数目，m为景观中的斑块类型总数目。

(8)河网密度，指1km×1km评价单元内河流的长度。提取步骤为：首先利用ArcGISPro对遥感解译的面状河流数据提取其中心线，再利用1km×1km的网格评价单元统计每个网格内的河流长度，最后转化成栅格数据，经统一的规范化处理后得到河网密度初始数据。

具体的，河流是人类生存、发展所需淡水资源的重要组成，是人类发展的资源基础和环境载体。河网的疏密能综合反映一个地区的自然地理条件，它常随气候、地质、地貌等条件不同而变化。

(9)人口密度，在分县人口统计数据的基础上，综合考虑了与人口密切相关的土地利用类型、夜间灯光亮度、居民点密度等多因素，将以行政区为基本统计单元的人口数据展布到空间格网上，从而实现人口的空间化，经统一的规范化处理后得到人口密度初始数据。

具体的，人口的激增对资源和环境带来了很大的压力。由于人口的增长，城乡不断扩展，大量的耕地被占用；而且过度利用土地导致耕地表土侵蚀严重，肥力下降；砍伐森林、开垦草原、围湖造田对生态平衡造成了严重的破坏。这些导致了人口增长与土地资源减少之间的矛盾越来越尖锐，人口增长对土地生态的破坏也越来越大。

(10)路网密度，指1km×1km评价单元内所有道路的总长度。提取步骤为：首先利用ArcGIS Pro对遥感解译的面状道路数据提取其中心线，再利用1km×1km的网格统计每个网格内的道路总长度，最后转化成栅格数据，经统一的规范化处理后得到道路密度初始数据。

具体的，道路是一种基础设施，在促进人类经济和社会发展的同时，干扰了原本处于稳定状态的生态系统，并产生了一系列的生态效应，直接或间接地影响到了生态系统的结构、功能和格局。

(11)不透水盖度，指1km×1km评价单元内不透水面的面积。提取步骤为：首先基于遥感解译的不透水数据，利用1km×1km的网格统计每个网格内不透水面面积，再转化成栅格数据，经统一的规范化处理后得到不透水盖度初始数据。

具体的，不透水面难以将降水渗透到土壤，导致地下水交换和基流活动减弱；而且不透水面具有较强的太阳辐射吸收能力，会使城市出现热岛效应。因此，不透水面盖度对流域的水热循环和水环境有着重要的影响。

(12)地均粮食产量，在粮食产量分县统计数据的基础上，提出一种基于MODIS GPP数据实现粮食产量空间化的模型，通过该模型计算得到地均粮食产量初始数据。

具体的，粮食生产是人类利用土地满足生存需求的社会经济活动，是人与自然和谐相处的结果。相关研究表明各地区作物产量与该地区作物生长季内的总初级生产力GPP关系密切，存在转换系数关系，因此利用MODIS GPP数据粮食产量统计数据进行空间化，详细提取方法如下：

首先，对全年GPP数据时间序列数据进行线性插值和S-G滤波，以及统一的规范化处理；

其次，提取每个耕地像元的全年GPP总和，记为GPP_T；

然后，利用县(市、区)域行政区矢量统计各县(市、区)范围内所有耕地像元的GPP_T的总和；

最后，利用像元GPP_T占县(市、区)域GPP_T的总和的比，再乘以县(市、区)域粮食产量，实现粮食产量的空间化，得到地均粮食产量初始数据。计算公式如下：

式中Y_ik为第k个县域行政区第i个像元的粮食产量，GPP_T_ik为第k个县域行政区第i个像元的GPP年总和值，u为第k个县域行政区的像元总数，Y_k为第k个县域行政区的粮食产量。

(13)地均GDP，在分县GDP统计数据的基础上，综合考虑与人类经济活动密切相关的土地利用类型、夜间灯光亮度、居民点密度等多因素，将以行政区为基本统计单元的GDP数据展布到栅格单元上，从而实现GDP的空间化，经统一的规范化处理后得到地均GDP初始数据。

具体的，推动中国经济高速增长的重要因素是土地要素。随着工业化和城市化的加速发展，土地价值在经济社会中越发受到重视，地方政府依靠土地优化资源配置，促进地方经济发展。

(14)人均生态用地面积，指1km×1km评价单元内生态用地面积与人口的比值。提取步骤为：首先基于遥感解译的生态用地，利用1km×1km网格统计每个格网内土地生态用面积，再转化成栅格数据，最后经统一的规范化处理后再除以人口密度数据，得到人均生态用地初始数据。

具体的，生态用地是指能直接或间接发挥调节气候、保持水土、净化空气、美化环境等生态功能并提供生物栖息地、维护生物多样性等生态保育功能，对促进区域生态平衡和区域持续发展具有重要作用的土地，主要包括耕地、林地、草地、水域和其他未利用土地等生态用地。

(15)生态服务价值，提取步骤为：首先计算1km×1km评价单元内不同地类面积和生态服务价值系数乘积的和，再经统一的规范化处理后得到生态服务价值初始数据，计算公式如下：

式中，ESV生态系统服务价值总量，单位：万元；A_s是研究区第s种土地利用类型分布面积；VC_s为第s种地类的生态服务价值系数。

具体的，以Costanza等的研究理论及谢高地等制定的价值系数为基础，得到不同地类的生态服务价值系数，见表2。

表2

土地利用类型	对应生态系统类型	生态服务价值系数/(元·hm-2)
			耕地	农田	3033.11
林地	森林	15278.72
			草地	草地	6154.57
水域	水体	33534.36
			城乡、工矿、居民用地	城市	88.57
未利用土地	荒漠	529.15

(16)植被覆盖度，是反映区域生态环境状况的指示因子。提取步骤为；首先对同年度区域所有遥感影像进行精确的辐射定标、大气校正和几何配准；然后利用云层指数CCI和阴影指数YYI消除云和阴影的影响，计算公式如下：

CCI＝(ρ_R+ρ_G+ρ_B)×ρ_NIR；

式中，ρ_R、ρ_G、ρ_B、ρ_NIR、ρ_SWIR分别为红波段、绿波段、蓝波段、近红外波段和短波红外波段的光谱反射率值，NDVI为归一化植被指数。

最后利用最大值合成法提取年度内所有NDVI数据的最大值，并采用像元二分法计算植被覆盖度，经统一的规范化处理后得到植被覆盖度初始数据。计算公式如下：

式中，VFC为像元的植被覆盖度，NDVI为该像元上的归一化植被指数，NDVI_soil为无植被覆盖像元的NDVI值，NDVI_veg为完全被植被覆盖的像元NDVI值。其中，NDVI_veg和NDVI_soil分别为NDVI累积百分数5％和95％对应的NDVI值。

(17)植被总初级生产力，首先对MODIS GPP时间序列数据进行线性插值和S-G滤波；然后计算年均值，经统一的规范化处理后得到研究区内的植被总初级生产力初始数据。

具体的，GPP是植物通过光合作用固定的碳总量，其时空变化趋势对土地以及全球碳源汇的优化计算至关重要。

S4，对步骤S3中的17个评价指标初始数据进行极差标准化处理，得到17个标准化后的评价指标极标数据；计算公式如下：

若评价指标为正向：

若评价指标为逆向：

式中，x_maxi、x_mini分别为第i个评价指标初始数据的最大值和最小值，x_mi为第i个评价指标初始数据在第m个评价单元的值，y_mi为第i个评价指标极标数据在第m个评价单元的值。

具体的，按指标对土地生态状况的影响特征不同，将各指标划分为正向指标、负向指标及区间指标。其中正向指标值越大地区土地生态状况越好，负向指标则相反，各指标性质见表1，其中地形位指数和年降水量常用作区间指标，以陡坡的最大值35°坡度值为阈值划分地形位指数生态效应临界值，由于本区域坡度均小于35°，所以将地形位指数作为正向指标；年降水量在以年份为时间尺度的情况下，难以判断其对土地生态状况的影响阈值或阈值区间，故也作为正向指标处理。

S5，针对步骤S4中评价指标极标数据，利用层次分析法和熵权法主客观综合赋权测算各评价指标的权重，见表3；

表3

S6，针对步骤S4中的评价指标极标数据和步骤S5中各评价指标极标数据的权重，利用综合指数法评估土地生态状况，得到河北省土地生态状况综合指数LESI，计算公式如下：

式中，LESI为土地生态状况综合指数；X_ij为第i个评价单元的第j个评价指标极标数据的值；w_j为评价指标极标数据X_ij的权重；z为评价指标极标数据的数量。

S7，针对步骤S6中的土地生态状况综合指数LESI，利用自然段点法划分土地生态状况等级，得到河北省土地生态状况等级空间分布结果；

其中，土地生态状况最好的区域主要位于太行-燕山一带。这是由于该区域生态基本状况、结构状况较好且受到的外界环境干扰少，分布有大量的具有重要生态功能的林地和草地，对应的植被覆盖度和GPP较优，生境调节能力强；同时在山地、高原地形因素和退耕还林、京津风沙治理等政策因素的共同影响下，该区域受人为干扰程度较低，土地利用结构、景观蔓延性较高，并带来了较高的土地生态服务价值。其次是坝下低中山盆地，该区域草地面积比例较高，植被覆盖度和GPP也相对较高，外界环境等人为干扰程度较低，并带来了较高的人均生态用地和土地生态服务价值。坝上高原、冀西北山区和平原区的土地生态状况普遍不高，其中坝上高原和冀西北山区的土壤条件、降雨量GPP等土地基础条件较差，虽然人类活动相对较少，但该地区土地结构稳定性差，土地利用带来的效益有限。而生态状况最低的区域主要分布于城市建成区及周边地区。这些区域的土地利用类型主要是耕地和建设用地，由于人类活动相对频繁，从而对土地生态状况的造成了较为显著的负面影响。

具体的，在ArcGIS中利用自然段点法将河北省土地生态状况划分差、较差、一般、较好、好五个等级。可以看出河北省土地生态状况形成了较为明显的地理梯度分布，土地生态状况好和较好的面积占比超过70％，但差和较差的土地状况区域面积占比超过了11％，土地生态状况仍有一定的提升空间。

尽管这里参照本发明的多个解释性实施例对本发明进行了描述，但是，应该理解，本领域技术人员可以设计出很多其他的修改和实施方式，这些修改和实施方式将落在本申请公开的原则范围和精神之内。更具体地说，在本申请公开、附图和权利要求的范围内，可以对主题组合布局的组成部件和/或布局进行多种变型和改进。除了对组成部件和/或布局进行的变形和改进外，对于本领域技术人员来说，其他的用途也将是明显的。

Claims (3)

1.一种耦合多源遥感信息的土地生态状况评估方法，其特征在于，包括以下步骤：

S1，基于17个评价指标构建多层级综合评价指标体系；

多层级综合评价指标体系包括生境状况、结构状况、胁迫状况和土地利用效益4个准则层，4个准则层包括10个要素层，10个要素层包括17个评价指标；

所述生境状况包括立地条件、土壤条件和水热条件；所述结构状况包括土地利用结构和生态景观；所述胁迫状况包括人口压力和人为干扰；所述土地利用效益包括社会效益、经济效益和生态效益；

其中，所述立地条件为地形位指数，所述土壤条件为土壤质地，所述水热条件包括年降水量和活动积温，所述土地利用结构为土地利用综合指数，所述生态景观包括景观多样性指数、景观蔓延度指数和河网密度，所述人口压力为人口密度，所述人为干扰包括路网密度和不透水盖度，所述社会效益为地均粮食产量，所述经济效益为地均GDP，所述生态效益包括人均生态用地、生态服务价值、植被覆盖度和植被总初级生产力；

17个评价指标分别为：地形位指数、土壤质地、年降水量、活动积温、土地利用综合指数、景观多样性指数、景观蔓延度指数、河网密度、人口密度、路网密度、不透水盖度、地均粮食产量、地均GDP、人均生态用地、生态服务价值、植被覆盖度、植被总初级生产力；

S2，采用1km×1km格网对研究区进行评价单元的划分；

S3，利用多源遥感数据并采用步骤S2中的1km×1km格网对步骤S1中的17个评价指标进行提取和优化，得到17个评价指标初始数据；

17个评价指标的提取和优化如下：

(1)地形位指数，利用规范化处理后的DEM遥感数据和坡度计算得到研究区内的地形位指数初始数据，计算公式如下：

式中，T为地形位指数；E和S分别表示区域范围内任意一点M的高程值和坡度值；和分别表示M点所在区域的平均高程值和平均坡度值，其中，坡度值由DEM数据通过ArcGISSlope工具计算得到；

(2)土壤质地，使用土壤粉砂含量，经统一的规范化处理后得到研究区内的土壤质地初始数据；

(3)年降水量，基于卫星遥感降水数据，对同年度区域所有日降水数据进行求和，经统一的规范化处理后得到研究区内的年降水量初始数据；

(4)活动积温，通过NDVI时间序列数据确定的植被生长开始日期和结束日期，计算植被生长开始日期至结束日期的地表温度LST累计和，经统一的规范化处理后得到研究区内的活动积温初始数据；

(5)土地利用综合指数，基于遥感解译的土地利用数据，从三生空间的角度出发，构建一种土地利用综合指数模型，并通过一种土地利用综合指数模型计算得到土地利用综合指数初始数据；

一种土地利用综合指数模型的构建包括以下步骤：

首先，从三生空间角度将遥感解译的土地利用数据分类为生态空间、生产空间和生活空间，其中生产空间是以提供工业品、农产品和服务产品为主导功能的区域；生活空间是以提供人类居住、消费、休闲和娱乐为主导功能的区域；生态空间是以提供生态产品和生态服务为主导功能的区域；

其次，对三生空间进行分级赋值，生态空间赋值3，生产空间赋值2，生活空间赋值1；

最后，利用1km×1km网格统计每个格网内各用地空间类型等级值和面积百分比乘积的和，并转化成栅格数据，得到土地利用综合指数初始数据，计算公式如下：

式中，L_c为土地利用综合指数，A_i为第i种用地空间类型的分级赋值，C_i为第i种用地空间类型的面积百分比，n为用地空间类型的总数量；

(6)景观多样性，使用Shannon Weaver指数度量，提取步骤为：利用最优尺度景观单元格网对遥感解译的土地利用数据进行测算，经统一的规范化处理得到景观多样性初始数据，计算公式如下：

式中，SHDI为景观多样性指数，p_i为景观斑块类型i所占的面积百分比，q为景观中的斑块类型总数目；

(7)景观蔓延度，提取步骤为：利用最优尺度景观单元格网对遥感解译的土地利用数据进行测算，经统一的规范化处理得到景观蔓延度初始数据，计算公式如下：

式中，CONTAG为景观蔓延度指数，p_i为景观斑块类型i所占的面积百分比，g_ik为i类型斑块和k类型斑块毗邻的数目，m为景观中的斑块类型总数目；

(8)河网密度，指1km×1km评价单元内河流的长度，提取步骤为：首先提取遥感解译的面状河流数据的中心线，然后利用1km×1km的网格评价单元统计每个网格内的河流长度，再转化成栅格数据，经统一的规范化处理后得到河网密度初始数据；

(9)人口密度，在分县人口统计数据的基础上，综合考虑了与人口密切相关的土地利用类型、夜间灯光亮度、居民点密度因素，将以行政区为基本统计单元的人口数据展布到空间格网上，从而实现人口的空间化，经统一的规范化处理后得到人口密度初始数据；

(10)路网密度，指1km×1km评价单元内所有道路的总长度；提取步骤为：首先提取遥感解译的县道、省道、国道和铁路数据的中心线，再利用1km×1km的网格统计每个网格内的道路总长度，最后转化成栅格数据，经统一的规范化处理后得到路网密度初始数据；

(11)不透水盖度，指1km×1km评价单元内不透水面的面积，提取步骤为：首先基于遥感解译的不透水数据，利用1km×1km的网格统计每个网格内不透水面面积，再转化成栅格数据，经统一的规范化处理后得到不透水盖度初始数据；

(12)地均粮食产量，在粮食产量分县统计数据的基础上，构建一种基于MODIS GPP数据实现粮食产量空间化的模型，并通过一种基于MODIS GPP数据实现粮食产量空间化的模型计算得到地均粮食产量初始数据；

一种基于MODIS GPP数据实现粮食产量空间化的模型的构建包括以下步骤：

首先，对全年GPP时间序列数据进行线性插值和S-G滤波，以及统一的规范化处理；

其次，提取每个耕地像元的全年GPP总和，记为GPP_T；

然后，利用县域行政区矢量统计各县范围内所有耕地像元的GPP_T的总和；

最后，利用像元GPP_T占县域GPP_T的总和的比，再乘以县域粮食产量，实现粮食产量的空间化，得到地均粮食产量初始数据，计算公式如下：

式中，Y_ik为第k个县域行政区第i个像元的粮食产量，GPP_T_ik为第k个县域行政区第i个像元的GPP年总和值，u为第k个县域行政区的像元总数，Y_k为第k个县域行政区粮食产量；

(13)地均GDP，在分县GDP统计数据的基础上，综合考虑与人类经济活动密切相关的土地利用类型、夜间灯光亮度、居民点密度因素，将以行政区为基本统计单元的GDP数据展布到栅格单元上，从而实现GDP的空间化，经统一的规范化处理后得到地均GDP初始数据；

(14)人均生态用地面积，指1km×1km评价单元内生态用地面积与人口的比值，提取步骤为：首先基于遥感解译的生态用地，利用1km×1km网格统计每个格网内土地生态用面积，再转化成栅格数据，最后经统一的规范化处理后再除以人口密度数据，得到人均生态用地初始数据；

(15)生态服务价值，提取步骤为：首先计算1km×1km评价单元内不同地类面积和生态服务价值系数乘积的和，再经统一的规范化处理后得到生态服务价值初始数据，计算公式如下：

式中，ESV为生态系统服务价值总量，单位：万元；A_s是研究区第s种地类分布面积；VC_s为第s种地类的生态服务价值系数；

(16)植被覆盖度，提取步骤为：首先对同年度区域所有遥感影像进行精确的辐射定标、大气校正和几何配准；然后利用云层指数CCI和阴影指数YYI消除云和阴影的影响，计算公式如下：

CCI＝(ρ_R+ρ_G+ρ_B)×ρ_NIR；

式中，ρ_R、ρ_G、ρ_B、ρ_NIR、ρ_SWIR分别为红波段、绿波段、蓝波段、近红外波段和短波红外波段的光谱反射率值，NDVI为归一化植被指数；

最后利用最大值合成法提取年度内所有NDVI数据的最大值，并采用像元二分法计算植被覆盖度，经统一的规范化处理后得到植被覆盖度初始数据，计算公式如下：

式中，VFC为像元的植被覆盖度，NDVI为该像元上的归一化植被指数，NDVI_soil为无植被覆盖像元的NDVI值，NDVI_veg为完全被植被覆盖的像元NDVI值，其中，NDVI_veg和NDVI_soil分别为NDVI累积百分数5％和95％对应的NDVI值；

(17)植被总初级生产力，提取步骤为：首先对MODIS GPP时间序列数据进行线性插值和S-G滤波；然后计算年均值，经统一的规范化处理后得到研究区内的植被总初级生产力初始数据；

其中统一的规范化处理为转投影、重采样和裁剪，以使17个评价指标初始数据的空间信息与1km×1km格网评价单元一致；

S4，对步骤S3中的评价指标初始数据进行极差标准化处理，得到17个标准化后的评价指标极标数据，极差标准化的计算公式如下：

若评价指标为正向：

若评价指标为逆向：

式中，x_maxi、x_mini分别为第i个评价指标初始数据的最大值和最小值，x_mi为第i个评价指标初始数据在第m个评价单元的值，y_mi为第i个评价指标极标数据在第m个评价单元的值；

S5，针对步骤S4中标准化后的评价指标极标数据，利用层次分析法和熵权法主客观综合赋权测算各评价指标极标数据的权重；

S6，针对步骤S4中标准化后的评价指标极标数据和步骤S5中各评价指标极标数据的权重，利用综合指数法评估土地生态状况，得到土地生态状况综合指数LESI，LESI计算公式如下：

式中，LESI为土地生态状况综合指数；X_ij为第i个评价单元的第j个评价指标极标数据的值；w_j为评价指标极标数据X_ij的权重；z为评价指标极标数据的数量；

S7，针对步骤S6中的土地生态状况综合指数LESI，利用自然段点法划分土地生态状况等级，得到土地生态状况等级空间分布结果。

2.如权利要求1所述的一种耦合多源遥感信息的土地生态状况评估方法，其特征在于，步骤S2中1km×1km格网评价单元是先以MODIS LST数据为环境变量，利用ArcGIS CreateFishnet工具生成格网；再利用研究区矢量进行裁剪得到。

3.如权利要求1所述的一种耦合多源遥感信息的土地生态状况评估方法，其特征在于，景观多样性和经蔓延度中的最优尺度景观单元格网确定方法如下：

首先，利用ArcGIS Create Fishnet创建1km、2km、2+N…50km共26个窗口尺度的景观单元网格；

其中，N为大于等于2且小于等于46的偶数；

其次，分别利用26个窗口尺度的景观单元网格裁剪土地利用数据，得到相应的景观单元的土地利用类型数据；

然后，将景观单元土地利用数据批量输入Fragstats4.2软件计算景观多样性指数SHDI；

最后，通过尺度效应曲线选取最优尺度的景观单元格网，用于提取景观多样性和景观经蔓延度。