融合时序遥感分析的国土空间生态保护修复关键区识别

陈元鹏 周 旭 陈 妍 刘岩涛 苏香燕 张成鹏

(1.自然资源部国土整治中心, 北京 100035; 2.中国地质大学(北京)土地科学技术学院, 北京 100083;3.南京大学地理与海洋科学学院, 南京 210023)

0 引言

国土空间生态保护修复是指遵循生态系统演替规律和内在机理,对生态功能退化、生态系统受损、空间格局失衡、自然资源开发利用不合理的生态、农业、城镇国土空间,统筹和科学开展山水林田湖草沙一体化保护修复的活动[1-3]。国土空间生态保护修复是落实国家生态文明建设的重要措施,也是构建国家生态安全格局和统筹山水林田湖草沙系统治理的重要举措[4-5]。作为提升生态系统结构和功能完整性的有效方案,国土空间生态保护修复已上升为国家战略工程[6]。

国土空间生态保护修复关键区识别至关重要,是国土空间生态保护修复规划编制、生态保护修复工程规划和时空布局、生态保护修复工程措施选取等系列工作的前置条件和基础保障,是国土空间生态保护修复工作的重点和难点[7-9]。

目前国土空间保护修复关键区的识别主要依赖于最小累积阻力模型[10-12]或电路理论[3, 13-14]构建的生态安全格局。“识别生态源地-构建阻力面-提取生态廊道”的框架模式成为当前生态安全格局构建的基本范式[15],其中生态源地的提取是重要基础[16]。生态源地是物种或生态流向外扩散的起点,可以促进生态过程、维持系统完整性和提供高质量的生态系统服务,是保障区域生态安全的关键区域[17]。识别生态源地的方法可以分为直接法[18]与综合评价法[19]。现在研究和采用较多的方法,是基于多指标生态系统服务并考虑人为因素识别生态源地。然而,目前的方案虽然充分考虑了变量间的数值关系,但对变量间的生态关系考虑不足,较少考虑不同生态系统服务之间的权衡协同场景[20-22]。

在目前生态安全格局构建的范式下,国土空间生态保护修复关键区识别的方法通常从单一时间节点出发,在识别生态源地的基础上提取生态廊道与生态战略点进而实现目标,但缺少对生态系统时间序列变化的考虑[23]。因此,本文在生态系统健康评估中考虑生态系统服务权衡协同的场景,更加科学合理地提取生态源地;另一方面在生态源地提取的基础上,增加对长时间序列植被指数变化趋势的分析,融合国土空间生态安全格局的动态和静态因素,提出一种更全面、准确、合理的生态保护修复关键区域的识别技术方法,为国土空间生态保护修复本底调查、问题识别、空间规划以及工程布局等提供技术支撑。

1 研究区与数据源

1.1 研究区概况

研究区为河南省某生态保护修复工程所在地,位于秦岭山脉东段北部、河南省西部,地理坐标为33°38′33″～35°04′50″N、110°21′50″～113°36′57″E。该区域地形地貌变化剧烈,由西向东形成秦岭中低山地、黄土丘陵台塬、伊洛山间盆地、黄河滩区湿地的过渡带地貌,第1级地貌为西部小秦岭、伏牛山等中山区,平均海拔1 000～2 000 m,部分山峰海拔超过2 000 m;第2级地貌为中部崤山、嵩山等低山丘陵一带,平均海拔200～1 000 m;第3级地貌为东北部伊洛盆地和黄河滩区,平均海拔200 m以下。研究区地形地貌如图1所示。

图1 研究区地形地貌图Fig.1 Topographic of study area

特有的区域气候和地貌特征奠定了实施区森林、湿地、河湖、农田、城镇等5类陆地生态系统发育与演变的自然基础,以及社会经济发展的空间格局。

研究区域主要为暖温带季风型气候,四季分明,年平均气温14.7℃,1月最冷,月平均气温0.3℃,7月最热,月平均气温27.5℃。年平均降水量606.9 mm,降水主要集中在6—9月,占全年降水量的62.4%。多年平均蒸发量为1 829.10 mm。

1.2 数据来源及处理

本文所用数据包括气象数据、遥感数据、土地覆被数据等。气象数据来源于中国科学院资源环境科学与数据中心网站,中国气象要素年度空间插值数据集,包括2020年气温、降水量、蒸发量等,空间分辨率为1 km。土壤数据来源于中国科学院资源环境科学与数据中心网站1∶1 000 000土壤数据库。地形数据为地面数字高程模型数据(DEM),来源于地理空间数据云SRTMDEM,空间分辨率为90 m。遥感数据为基于MOD13A1构建的时序(2011—2020年)归一化植被指数(NDVI)数据集,共计230景,空间分辨率为500 m。土地覆被数据为MODIS数据产品MCD12Q1土地覆被类型(Land cover),选择其中2010、2015、2020年3期,空间分辨率为 500 m。 植被净初级生产力(NPP)来源于MODIS数据产品MOD17A3HGF数据集,时间序列为2011—2020年,空间分辨率为500 m。其中,时序NDVI、NPP、土地覆被类型数据均通过Google Earth Engine(GEE)云平台获取并批量处理计算。

GEE是当前世界上较为先进的处理卫星影像等地理空间观测数据的云计算平台。GEE云端数据库中集成了近40年的Landsat系列卫星的历史存档数据,给个人用户提供了强大的算力和云存储空间,同时提供了方便快捷的JavaScript语言API接口进行数据处理、算法实现以及结果分析[24-25]。本研究中,通过GEE云计算平台对研究所需的遥感、土地覆被、NPP数据进行了批量处理,减少了前期大量数据的准备工作,有效提升了工作效率,降低了数据处理与算法实现过程中对本地硬件设备的依赖。

2 研究方法

本文研究方法框架主要包括3部分:首先,基于GEE云平台、MOD13A1数据,采用Mann-Kendall方法对研究区2011—2020年时序NDVI开展趋势分析,形成时序分析结果;其次,基于气象、地形、土壤、NPP等数据,开展基于生态系统服务协同和权衡计算的生态系统健康分析评估,选择生态源地;最后,基于时序分析结果和生态源地选择结果,采用叠置分析,开展生态保护修复关键区域识别。具体技术路线如图2所示。

图2 技术流程图Fig.2 Workflow of this study

2.1 Mann-Kendall趋势分析(时序遥感趋势分析)

Theil-Sen Median(TSM)趋势分析[26]和Mann-Kendall(MK)非参数检验法[27]常用来分析NDVI空间演变特征。本文采用 Theil-Sen Median 趋势分析方法并结合MK检验方法判断Sen趋势的显著性,与一元线性回归趋势分析方法相比,其可以规避时间序列数据分布和数据缺失对分析结果的影响[28],同时可以降低异常值对结果影响[29]。TSM计算公式为

(1)

式中i、j——时间序列数

NDVIi、NDVIj——第i、j时间序列的NDVI值

Median——多年NDVI中位数函数

n——时间序列长度

MK是一种非参数统计检验方法,用来判断趋势的显著性,其计算公式为

(2)

(3)

本文使用显著性水平α=0.05,对应的u1-α/2为1.96[30]。

2.2 基于生态系统健康的生态源地识别

本文在生态系统服务权衡与协同框架下开展生态系统健康评估,进而开展生态源地识别。

2.2.1生态系统健康评估体系

生态源地是对区域生态过程和功能起决定性作用的生境斑块,提供必要的生态系统服务,对区域生态系统健康安全具有重要意义。生态系统健康(Ecosystem health,EH)评估是生态源地识别的基础[31],本文从生态系统活力(Ecosystem vigor)、生态系统结构(Ecosystem organization)、生态系统弹性(Ecosystem resilience)、生态系统服务(Ecosystem service,ES)4个维度构建评估框架体系[32]。计算公式分别为

(4)

(5)

式中EH——生态系统健康水平

PH——生态系统物理健康水平

ES——生态系统服务

V——生态系统活力

O——生态系统结构

R——生态系统弹性

ES、V、O和R计算结果均进行归一化处理。

V可用NPP表征,因为NPP能够很好地反映植物固定和转化光合产物的效率,并表征可供人类使用的物质和能量。

O表示区域生态系统结构的稳定性,主要包括景观异质性和连通性。为进一步描述景观异质性和连通性,选择景观格局指数(包括面积加权平均斑块分维数(AWMPFD)、香农多样性指数(SHDI)、景观聚集度指数(CONT)、景观破碎化指数(FN))用以表征O,公式为

O=0.3SHDI+0.2AWMPFD+
0.25FN+0.25CONT

(6)

式(6)中,景观格局指数采用Fragstats软件计算。

R表示生态系统在面对人为扰动时保持自身结构稳定性的能力,反映一个区域在生态系统过程中抵御和适应外部干扰的能力,因此包括两方面,即抵抗力和恢复力[33],本研究中为抵抗力和恢复力分别赋值0.6和0.4[34],为反映同一土地覆被类型不同生态系统的弹性,采用NDVI进行校正,公式为

(7)

式中Ri——第i个像元总弹性系数

NDVIm——第m个像元的NDVI值

NDVImeanj——第j类土地覆被类型的NDVI均值

Resilj、Resisj——第j类土地覆被类型的恢复力系数和抵抗力系数(表1)

表1 生态系统弹性系数Tab.1 Ecological resilience coefficient of land use types

生态系统服务是自然系统和社会经济系统之间的桥梁[35],反映了人与自然之间的耦合关系,是生态系统健康的重要组成部分。针对研究区的自然特征和生态环境状况,本文选取水源涵养、水土保持、防风固沙3类生态系统服务作为区域关键生态系统服务功能进行定量评估。3类生态系统服务按照净初级生产力(NPP)定量指标评估法计算,计算公式分别为

WR=NPPmeanFsicFpre(1-Fslo)

(8)

Spro=NPPmean(1-K)(1-Fslo)

(9)

Sws=NPPmeanKFqD

(10)

其中

式中WR——生态系统水源涵养服务能力指数

NPPmean——多年植被净初级生产力平均值

Fsic——土壤渗流因子

Fpre——多年平均降水量因子

Fslo——坡度因子

Spro——生态系统水土保持服务能力指数

K——土壤可蚀性因子

Sws——生态系统防风固沙服务能力指数

Fq——多年平均气候侵蚀力

u——高度2 m的月平均风速

ETPi——月潜在蒸发量

Pi——月降水量

d——当月天数

Ti——月平均气温

ri——月平均相对湿度

D——地表粗糙度因子

θ——坡度

2.2.2生态系统服务权衡与协同

在生态系统服务的背景下,生态源地的选择过程必须充分考虑生态系统服务的权衡和协同关系。目前,有序加权平均模型(Ordered weighted average model,OWA)已有效用于量化生态系统服务之间的权衡和协同关系[36],因此基于OWA可量化不同场景下的生态系统服务供给,反映不同生态系统服务之间的权衡和协同效应[34, 37],其公式为

(11)

式中axj——第x栅格数据中第j个像元的属性值

Sxj——重新赋权后的栅格数据

ωx——Sxj的权重

N——参与计算的生态系统服务栅格数据数量

此外,采用风险值(risk)和权衡值(trade-off)两个参数作为OWA的约束和优化条件,公式为

(12)

(13)

(14)

通过计算风险值和权衡值,可以得出不同场景模式下的最优化的生态系统服务权重系数组合。计算最大权衡值和相应的权重,共设定11个场景。

2.2.3生态源地识别

为有效集约地选择生态源地,对于不同的场景模式,选择生态系统健康值前60%的区域作为优先保护区域,并采用生态系统保护效率公式在优先保护区域中进一步筛选生态源地。公式为[38]

(15)

2.3 生态保护修复关键区识别

以时序遥感趋势分析中NDVI的下降趋势、保持稳定和轻微上升趋势,以及生态系统健康评估和生态源地识别结果为基础,采用叠置分析,进一步识别生态保护修复关键区[39]。

3 结果与分析

3.1 时序遥感趋势

基于GEE云平台得到了研究区域2011—2021年的NDVI时间序列数据,共计230景,从统计描述看,230景NDVI数据的均值在0.16～0.76范围内波动,波动周期符合年度季相特征。标准差在0.06～0.23范围内波动,稳定性较好,数据质量较好,时序NDVI数据构建结果符合预期。研究区数据统计见图3。

图3 时序NDVI数据的均值和方差Fig.3 Mean and variance of time-series NDVI data

基于TSM趋势分析,结合Mann-Kendall非参数检验法,在显著性水平α=0.05下,形成时序NDVI趋势变化统计数据和空间分布(图4)。首先,从NDVI趋势变化统计数据结果看,研究区域内时序NDVI的斜率变化范围为-0.001 7～0.002 1,均值为0.000 3,其中斜率为负值的像元6 136个,占比5.33%,说明2011—2020年间,研究区域内少数局部区域NDVI呈下降趋势,但整体NDVI变化呈上升趋势,一定程度反映出研究区域内的生态系统质量稳中向好。其次,叠加土地覆被类型(图5)看NDVI趋势变化空间分布可知,2011—2020年间研究区域内NDVI呈下降趋势的区域主要分布在中部洛河、伊河沿线的不透水面集中地区及周边缓冲区域,东北侧黄河沿线的不透水面集中地区及周边缓冲区域,以及西北角小秦岭北侧的山前平原,其余区域还有些零星分布。NDVI呈上升趋势主要分布在中部黄土丘陵和东部山间盆地区域,主要土地覆被类型为耕地、草地和少量灌木林地等。NDVI变化趋势不显著的区域主要分布在秦岭中低山地区域,包括崤山、熊耳山、伏牛山等山脉,其土地覆被类型以林地为主。

图5 2020年土地覆被类型Fig.5 Spatial distribution of land cover types in 2020

3.2 生态系统健康评估结果

3.2.1生态系统物理健康(PH)

根据式(5)～(7),分别计算生态系统活力(V)、生态系统结构(O)、生态系统弹性(R),并经归一化处理,进一步计算得到生态系统物理健康(PH)评估结果及其空间分布(图6)。如图6所示,生态系统物理健康均值为0.592,研究区域内高值区主要分布在西侧秦岭中低山地区域,如小秦岭、熊耳山、伏牛山等,中低山地区域植被总体生长良好,生态系统健康程度高。由西向东,生态系统物理健康值逐渐下降,低值区域主要分布在伊洛山间盆地、河流滩区等土壤侵蚀严重、水土流失严重的区域。此外,城市化程度较高的地区,由于人类活动的干扰较大,生态系统物理健康水平较低。

图6 生态系统物理健康(PH)评估结果和空间分布Fig.6 Evaluation result and spatial distribution of ecosystem physical health

3.2.2生态系统服务(ES)

根据式(8)～(10)净初级生产力(NPP)定量指标评估算法,并经归一化处理,计算得到生态系统服务(ES)评估结果及其空间分布(图7)。不同生态系统服务能力之间存在空间异质性,反映了生态系统服务能力之间的权衡性,也揭示了不同生态过程对于研究区整体生态安全格局的影响和作用。如图7所示,研究区域内水源涵养服务能力由西南向东北逐步减弱,水源涵养能力较强的区域主要分布在中低山区,但并非区域内所有中低山区都具有较强的水源涵养服务能力,较强区域主要集中在伏牛山-熊耳山-外方山区域。较于水源涵养能力,水土保持能力分布相对均衡,重要区域主要分布在中低山、黄土丘陵区以及陆浑水库周边区域。防风固沙能力与水源涵养能力存在一定协同性,较强区域主要集中在中低山区,其中最强的区域主要分布在小秦岭、熊耳山、伏牛山区域,次强区域为崤山、青要山、外方山和嵩山区域。从单一生态系统服务空间格局看,研究区域内生态系统服务能力较强的区域主要分布在中低山区、黄土丘陵区等海拔较高且植被覆盖较强的区域。

图7 生态系统服务(ES)评估结果和空间分布Fig.7 Evaluation result and spatial distribution of different ecosystem services

3.2.3生态系统服务协同与权衡

根据式(11)～(14)生态系统服务权衡与协同算法,通过改变决策风险和协同权衡值,可以生成一系列场景。结合工作实际,以0为起始、0.1为间隔设置风险值,计算最大协同权衡值及其相应的权重,经计算共得到11个场景(表2)。

表2 不同风险场景及权衡水平下的最优系数Tab.2 Optimal order weights for different decision risk scenarios and trade-off levels

3.2.4生态系统健康(EH)评估和生态源地识别结果

根据生态系统物理健康(PH)和生态系统服务(ES)计算结果,结合生态系统服务权衡与协同的11个场景,按照式(4)计算得出11个生态系统健康(EH)评估结果及其空间分布特征(图8)。

图8 生态系统健康(EH)评估结果及其空间分布Fig.8 Evaluation result and spatial distribution of ecosystem health in different scenarios

如图8所示,生态系统健康(EH)的高值和低值空间分布呈一致趋势,但局部的高值区和低值区也存在一定差异,这是由于不同的场景、不同的生态系统服务权衡与协同水平下,权重等次不同,一定程度上说明了简单的加权求和不能充分反映区域生态系统健康(EH)的空间分布,因此,有必要将生态系统服务协同权衡效应纳入生态系统健康评估。

基于11个场景的生态系统健康评估结果,按照高值的前60%筛选了不同场景下的优先保护区域,并根据式(15)测量了其保护效率(表3)。如表3所示,场景9的保护效率最高,为3.02,因此选择场景9的高值前60%区域作为生态源地,其空间分布如 图9 所示。研究区域内生态源地约1 386.11 km2,主要分布在西南部的伏牛山-熊耳山-外方山区域,少量分布在小秦岭区域,涉及的土地覆被类型主要为林地和草地。

表3 不同场景下的优先保护区域保护效率EiTab.3 Protection efficiency of priority protected areas in different scenarios

图9 生态源地选择结果及其空间分布Fig.9 Select result and spatial distribution of ecological sources

3.3 生态保护修复关键区

以时序遥感趋势分析结果(图4)为基础,选取NDVI下降趋势(斜率小于0)、NDVI保持稳定和轻微上升趋势(斜率在0～0.000 2之间)的区域,与生态源地选择结果(图9)进行叠置分析。

将NDVI保持稳定和轻微上升的区域和生态源地合并为生态保护关键区,用以表征研究区域内需生态保护的关键区域,其中:①将NDVI保持稳定和轻微上升区与生态源地重合的区域划定为一级生态保护区,因为该区域既是保持生态系统健康的重要区域,但从时序上看,近10年间该区域可以表征生态系统的NDVI指数并未形成显著变化趋势,一直属临界状态,如不作为重点保护区域则可能出现生态系统退化或破坏情况,因此该区域有必要作为一级生态保护区,以逐步形成生态系统稳定性较强、NDVI指数稳中有升的区域。②将NDVI上升区与生态源地重合的区域划定为二级生态保护区,二级生态保护区同样是保持生态系统健康的重要区域,但与一级生态保护区相比,近10年间区域内NDVI已形成上升趋势,生态系统已稳中向好,故此作为二级生态保护区。③将未与生态源地重合的NDVI保持稳定和轻微上升区划定为三级生态保护区,与一、二级生态保护区相比较,三级生态保护区未包含生态源地,其对于保持区域内生态系统完整、提供高质量生态系统服务的贡献度相对不强,但近10年间该区域内生态系统未遭到明显破坏,因此有必要进行保护以维持现状或将其逐步转化成为生态系统稳定性较强、NDVI指数稳中有升的区域。④将NDVI下降区划定为生态修复区,采取适当的措施实施生态修复。形成的生态保护修复关键区识别结果及其空间分布如图10所示。

图10 生态保护修复关键区识别结果及其空间分布Fig.10 Spatial distribution of ecological protection and restoration key areas identification results

由图10可知,一、二级生态保护区主要分布在西部和西南部小秦岭、熊耳山和伏牛山等区域,主要为中山森林生态系统,承担了生物多样性保护、固碳释氧等生态系统服务功能;三级生态保护区主要分布在邙山、崤山区域,在伊洛盆地和黄河滩区也有零星分布,主要为低山丘陵森林生态系统、山间盆地农田生态系统、河流生态系统,承担了水源涵养、水土保持、固碳释氧、矿产资源提供、农林产品提供生态系统服务功能;生态修复区主要分布在中部洛河、伊河沿线和下游不透水面集中地区及周边缓冲区域,东北侧黄河沿线的不透水面集中地区及周边缓冲区域,以及西北角小秦岭北侧的山前平原,主要为盆地区城镇、农田、河流生态系统,承担了农产品提供、城镇生态涵养生态系统服务功能。

4 结束语

基于时序遥感分析、生态系统服务协同权衡计算、生态系统健康评估、生态源地识别和空间叠置分析等方法,对研究区的生态保护修复关键区进行了有效识别。结果表明,提出的“基于生态系统服务协同权衡计算的生态系统健康评估——时序遥感趋势分析”的研究框架,对于识别构建区域生态保护修复关键区尤为重要。结合生态系统服务协同权衡的生态系统健康评估能够更有效地选择关键生态源地,在此基础上创新性地结合时序遥感分析进一步识别生态保护修复关键区,不仅顾及了研究区域内“静态”的生态系统服务、生态系统健康属性,同时衡量了“动态”的生态系统变化趋势,更全面、准确地识别了研究区的生态保护修复关键区域。研究成果可为国土空间生态保护修复本底调查、问题识别、规划和工程布局等提供技术支撑。