基于土地利用类型转移的甘肃祁连山自然保护区土地生态安全研究*

闫新宇, 牛成英, 任潇潇

(兰州财经大学统计学院,730020,甘肃省兰州市)

0 引 言

祁连山是我国西部重要的生态安全屏障,是黄河流域重要的水源产流地,也是我国生物多样性保护优先区域,加强祁连山生态系统整体保护对筑牢生态安全屏障具有重要意义. 2019年习近平总书记甘肃考察期间强调了祁连山生态保护的作用与生态屏障的作用,认为祁连山生态保护是国家的战略格局[1]. 由此可见对保护区生态安全及时作出评价,分析影响生态安全因素,对生态系统整体保护和系统修复具有重要指导意义.

土地生态安全是评价一个国家和地区生态安全的重要组成部分,然而,土地利用数据不明晰,林地、草地面积空间重叠,各部门统计数据不统一等问题是土地生态安全研究面临的一个重大困难[2]. 因此,基于遥感影像数据资料,对土地利用时空变化进行分析,从而研究土地生态安全影响因素是研究者常用的方法之一[3,4].

已有的文献对祁连山土地利用研究多为从空间地理或生态学方法上进行分析,主要研究土地利用对生态系统功能价值和地区气候环境的影响,或针对某一地类进行分析,而对土地利用与土地生态安全关系的研究尚显不足. 近年来,研究者越来越关注影响生态安全的多维度因素,通过分析影响生态安全的因素,构建生态安全指数(ESI)[5],从而对生态安全进行定量测度研究[6,7]. 目前,该指数已被广泛用于生态安全评价的多个领域[8,9],土地生态安全评价中生态安全指数也是常用方法之一[10,11],并且为了研究影响生态安全的影响因素,提出利用灰色关联分析方法分析生态安全关联因素[12,13].

而在面对资源环境发展的困境下,以环境保护为目的的各级环境资源规制相继出台,使得基于环境规制视角下的生态安全评价与保护引起研究者的关注[14]. 与此同时,国土资源规制在我国的实践与完善在各地试点展开,包括法律、行政、经济与技术等手段,如何借助国土资源手段对甘肃省祁连山自然保护区土地生态安全进行定量客观分析值得我们探索研究. 本文借助国土资源规制的技术手段,基于甘肃省祁连山自然保护区1985～2019年的土地卫星遥感数据,利用ArcGIS软件分析了祁连山自然保护区土地利用的时空转移特征,用文献研究分析法构建了基于土地利用类型的土地生态安全指数,利用灰色关联分析方法,对不同类型土地转移量与生态安全指数进行灰色关联分析,分析影响保护区土地生态安全的土地利用类型变化特征,为祁连山自然保护区的土地生态修复工作提供量化分析参考.

1 研究区域与数据来源

1.1 研究区域

甘肃祁连山自然保护区地处青藏、蒙新、黄土3大高原交汇地带的祁连山北麓,位于甘肃省境内祁连山北坡中、东段,地跨武威、金昌、张掖3市的凉州、天祝藏族自治县、古浪、永昌、甘州、山丹、民乐、肃南裕固族自治县8县(区),总面积为265.3万hm2,属森林生态系统类型的自然保护区. 本文研究的区域主要包括祁连山自然保护区的核心区、缓冲区和试验区,地理位置位于北纬36°93′～39°61′,东经97°40′～103°92′,如图1所示.

图1 甘肃省祁连山自然保护区分布图

1.2 数据来源与处理

本文的研究数据主要来源于国家青藏高原科学数据中心(http://data.tpdc.ac.cn/zh～hans/)和中国科学院资源环境科学数据中心(http://www.resdc.cn)数据库. 基于以上数据库中祁连山自然保护区30 m×30 m区域土地覆盖分类产品遥感数据集,选取1985～2019年间共11期土地利用数据,分别为1985、1990、1995、2000、2005、2010、2015、2016、2017、2018、2019年土地利用数据. 对祁连山自然保护区土地利用类型时空转移特征进行分析,进而分析土地利用变化对土地生态安全的影响因素.

首先,在数据提取与处理过程中,上述数据库中的数据集基于Google Earth Engine平台的Landsat系列数据经过分类处理后,获取了1985～2015年每5年一期的土地覆盖分类产品,以及2016～2019年每1年期的土地利用分类产品. 数据集的地理坐标系均为D_WGS_1984.

其次,对于1985～2019年祁连山土地利用栅格数据构建栅格属性表,依次对祁连山自然保护区的矢量数据按照掩膜提取,得到了1985～2019年祁连山自然保护区土地利用栅格数据. 基于保护区栅格数据的地理坐标系进行各土地利用类型面积的测算,得到1985～2019年祁连山自然保护区9类土地利用类型面积数据,包括未分类土地、耕地、林地、灌丛、湿地、水体、建筑用地、裸地、冰川/积雪9类. 根据常见分类,将上述9类土地利用类型整合为6类：耕地、林地(林地和灌丛)、草地、水体(湿地、水体、冰川和积雪)、建筑用地、未使用土地(裸地).

2 研究方法

2.1 土地利用转移矩阵及其基本原理

土地利用转移矩阵是对二维联列表的改造. 它将联列表的行列2个变量设置为2个时点,描述2个时点内的各土地利用类型间的相互转化. 由于在一个封闭的系统内,T1至T2时间段,各土地利用类型的转化是一个相互的过程,即系统内一种土地利用类型的增加必然会导致另一种土地利用类型的减少. 土地利用转移矩阵如表1所示[15].

表1 土地利用转移矩阵

其中,Pij表示T1至T2期间土地利用类型i转换为土地利用类型j的量占土地总量的百分比；Pi·-Pjj为T1至T2期间土地利用类型i转出为其他土地利用类型占土地总量的百分比,也就是T1至T2期间变量i减少的百分比,称为转出(Loss)；P·j-Pjj为T1至T2期间其他土地利用类型转化为土地利用类型j的量占土地总量的百分比,定义为T1至T2期间变量j增加的百分比,称为转入(Gain). 净变化量Dj为

Dj=max(Pj·-Pjj,P·j-Pjj)-min(Pj·-Pjj,P·j-Pjj)=|P·j-Pj·|.

(1)

在土地利用转移分析中,净变化量反映的是土地利用类型数量上的绝对变动,可以直观感受各土地利用类型面积的变动. 但土地利用是空间上不同地块间的转化,因而仅考虑净变化并不能真实地反映土地利用的动态演变.

为了避免此类情况,在净变化量的基础上引入了土地利用的交换变化量. 第j类土地利用类型交换变化量为

Sj=2min(Pj·-Pjj,P·j-Pjj),

(2)

而第j类土地利用类型的总变化量为

Cj=Dj+Sj=max(Pj·-Pjj,P·j-Pjj)+min(Pj·-Pjj,P·j-Pjj).

(3)

2.2 灰色关联分析

灰色关联度分析是一种多因素统计分析方法. 通过灰色关联度来描述因素之间发展趋势的相似或相异程度,进而衡量因素间的关联度[16].

假设有k个序列,每个序列有n期,灰色系统关联分析的具体计算步骤如下：

(1)根据方案评价的目的选取参照序列,记作

X0(t)=(x0(1),x0(2),…,x0(n))，

关联分析中被比较的序列常记为

Xi(t)=(xi(1),xi(2),…,xi(n)),i=1,2,…,k.

(2)数据标准化处理,正负向指数分别为

(3)计算关联系数

(4)

其中,Δi(t)=|X0(t)-Xi(t)|,i=1,2,…,k,t=1,2,…,n,Δi(t)min=min|X0(t)-Xi(t)|,Δi(t)max=max|X0(t)-Xi(t)|,ρ为分辨系数,取值范围为(0,1). 其关联度为

(5)

将各土地利用类型与参照序列生态安全指数之间的灰色关联度作为评价影响程度的度量值,通常认为灰色关联度越大的对象其影响程度越强.

3 祁连山自然保护区土地利用基本特征及生态安全指数变化趋势

3.1 土地利用基本特征

根据上述数据处理方法,可得甘肃省祁连山自然保区的土地利用类型面积及占比. 表2为耕地、林地(林地和灌丛)、草地、水体(湿地、水体、冰川和积雪)、建筑用地、未使用土地(裸地)6类土地利用类型占地面积与占比.

表2 1985～2019年土地利用面积(km2)与占比(%)

由表2可以看出,1985～2019年期间,祁连山自然保护区土地利用类型中草地占比最多,占比超过总面积50%以上,其次是未使用地,占比超过20%以上,占比最少为建筑面积,约为总面积的0.01%. 1985～2019年,土地利用变化幅度最大的为水体面积与未使用土地面积,耕地与林地面积变化基本处于稳定状态,耕地占比在区域总面积的2.4%～2.6%之间,林地占比在区域总面积的7.3%～7.5%之间. 1985～2005年期间,草地占比基本稳定在区域总面积的52%～53%之间. 2000～2010年间,草地变化较大,期间覆盖面积增加了约2%,其后在2010～2019年间存在小幅波动,占比基本稳定.

利用ArcGIS软件分析,1985～2019年间,祁连山自然保护区6类土地利用面积空间变化状况,可以直观的看出,祁连山自然保护区大体呈现为西北自东南的狭长走向,草地和裸地面积是研究区内的主要土地利用类型. 从方位上看,西北方向上的土地利用类型多为山体裸地和高原冰川,东南方向纬度相对较低,存在着广阔的草地和部分耕地,林地则多分布于研究区域的中心地区. 对比1985、1995、2005、2015与2019年各土地利用类型的空间分布情况可见,水体与未使用土地面积在时间和空间上均存在显著变化.

3.2 土地生态安全指数及其变化趋势

生态安全指数是伴随着生态学与地理学对于指数的研究发展而来,具体由生态环境指数演化而来.为探讨土地利用面积变化对区域土地生态环境的影响,根据实际情况,参考相关文献,测算祁连山自然保护区1985～2019年每一期的土地生态安全指数(ESI)[8,9]

(6)

其中,n为土地利用类型数,Ai为第i类土地利用面积,A为研究区总面积,Wi为第i类土地利用生态风险强度系数.

通过测算某一时间段内土地生态安全指数,可以综合评价研究区域土地生态安全程度及其变化趋势. 为了探究土地利用类型变动特征对土地生态安全的影响,本文分析近20多年祁连山自然保护区土地生态安全变化特征,采用文献研究分析法确定土地利用的生态风险权重[17-19],取不同土地利用类型生态风险系数的平均值为本文研究土地生态风险强度系数.

表3 不同土地利用类型的生态风险强度系数

由式(6)测算得祁连山自然保护区1985～2019年间11期土地利用生态安全指数,如表4所示.

表4 1985～2019年祁连山自然保护区土地生态安全指数

由表4结果可见,1985～2019年,祁连山自然保护区土地生态安全指数始终保持在0.68～0.72之间,土地生态安全指数平均水平处于较安全水平. 但由此也可以看出,1985～2019年,祁连山自然保护区土地生态安全出现了较大变化. 20多年来,祁连山自然保护区土地生态安全指数呈先下降,后上升的波动趋势. 1985～2010年间,祁连山自然保护区土地生态安全指数持续下降,说明这一时期土地生态安全持续衰退；2010～2016年间土地生态安全指数变化不显著；2016年之后,土地生态安全指数表现为在轻微波动中持续向好的态势,反映了该时期土地生态安全正在逐步恢复.

4 土地利用类型时空转移特征与生态安全指数灰色关联分析

4.1 土地利用类型时空转移特征

根据上述土地生态安全指数变化趋势,可将1985～2019年生态安全划分为3个状态,分别为1985～2010年的衰退期、2010～2015年间的涵养期、2015～2019年的修复期.

4.1.1 衰退期(1985～2010年)

这一时期生态安全指数持续降低,由1985年的0.7131降低到到2010年的0.6908. 由式(1)～式(3)计算得到土地利用转移净变化量,对这一时期土地利用时空转移特征进行分析可知,衰退期水体面积的变化最为突出,总变化量占总体变动的47.946%,具体表现为：水体面积转出9.895%,转入0.156%,净转出9.739%；其次是裸地面积变化较大,占总体变动面积的31.585%,转入6.135%,转出0.487%,净转入5.648%. 从以上分析可以看出,衰退期保护区内水体面积减小,裸地面积增加.

4.1.2 涵养期(2010～2015年)

2010～2015年期间,土地生态安全指数由2010年的0.6908变动到2015年的0.6887,变化幅度很小,可以定义为涵养期. 这一时期土地利用转移净变化量看出,涵养期的土地利用总变化量为5.234%,不足衰退期总变动率的1/4,涵养期6种土地利用类型的变动趋势基本稳定. 这一时期的裸地和水体面积分别净转出了0.886%和0.053%,同时草地面积净转入0.786%.

4.1.3 修复期(2015～2019年)

自2015年之后,祁连山自然保护区的土地安全指数总体表现出了逐步升高的趋势,由2015年的0.6887上升到了2019年的0.6967,因而,这一时期可以定义为土地修复期. 由这一时期的土地利用净变化量可见：修复期依旧是裸地与水体的变化最为突出,分别占总变动的45.238%和47.728%,同时草地也有较为明显的变动,为总变动的6.703%.

从用地类型转入的角度来看,修复期水体的转入主要来自于裸地和草地,其比重分别占总地类的5.054%和0.609%,与衰退期的表现形成了极大的反差. 同时也有0.017%的耕地退耕还林转化为草地.

综合上述3个变化时期土地利用类型时空转移特征,绘制地理空间变化图形(见图2)直观显示祁连山自然保护区土地利用时空演变状态.

图2 3个时期祁连山自然保护区土地利用时空演变图

由图2可直观看出,与1985年相比,2019年祁连山自然保护区的水体面积变化显著,原部分冰川转化为裸地和草地的痕迹较为明显. 由此可见,水体面积降低的最要因素是冰川/积雪面积的缩小. 1985～2019年,祁连山脉附近出现了由高原冰川到裸地再到草地或水体的转变,当下高原冰川的总面积仍然低于1985年的水平. 同时在研究区东南方向上的低纬度地区出现了由裸地到草地再到林地的良性转变,研究区林地面积相较于1985年出现了0.033%的增长.

4.2 各类土地利用类型转出与生态安全指数灰色关联分析

为了研究土地利用类型变化对生态安全指数波动的影响原因,本文基于ArcGIS软件,使用相交分析工具对上述6类共11期的土地利用数据进行了土地转移分析,得到了10个时期间6类土地利用类型的转出值,如表5所示. 由于祁连山自然保护区地理特征,表中各期的建设用地转出均接近于0,这说明期间建设用地变动相对较小. 由表5可见,1985～2019年间土地利用类型变动的主要体现在裸地和水体的转出.

表5 逐期起点与终点相比6类土地利用类型转出面积 单位：km2

为研究土地利用类型时空演变对土地生态安全的影响,有必要关注各类土地利用类型转出面积与土地生态安全指数的关联关系. 本文以土地生态安全指数为基期参照序列,以裸地、耕地、林地、草地、水体为比较的序列,讨论土地利用类型转出量与土地生态安全指数的关联水平. 其中,考虑到裸地和耕地的增加明显会对土地安全存在负向影响,裸地和耕地面积转出越多越有利于提高生态安全水平,因此将裸地、耕地转出量记为正向指标；其余土地利用类型转出量记为负向指标,即表示转出越多,越不利于地区生态安全.

根据灰色关联分析方法基本思想,有式(4)～式(5)可以得5类土地利用类型的转出量与土地生态安全指数的灰色关联度,如表6所示.

表6 各土地利用类型转出量与生态安全指数的灰色关联度

由表6可以看出,草地与水体转出量与土地生态安全指数关联度较高,分别为0.7343和0.7269,说明这2种土地利用类型的转移变化对土地生态安全指数存在较强的关联性,因而草地与水体面积的转出是导致土地生态安全指数波动的主要原因. 但宏观来看,草地的转出量较小,且呈现稳步的上升趋势,而水体转出量最大,且呈现不稳定波动趋势；耕地与土地生态安全指数也有较高的关联性,达到0.6428,但从本次研究的数据库提供数据来看,整个研究时间段林地的转出面积很小,因此,尽管与生态安全指数的关联强度达0.6591,但对土地生态安全变化的影响较小；裸地关联强调较低,但长期看,裸地的转出对生态安全存在正向影响,即裸地转出越多,生态安全指数越高.

5 研究结果分析

从上述土地利用类型转移特征及其与土地生态安全指数的关联分析结果可见,土地利用类型的转移对祁连山自然保护区土地生态安全指数的影响作用较为复杂,在每一时间段生态安全指数的变化受几种土地利用类型变化的综合作用.

(1)1985～2010年的衰退期. 从用土地利用类型转出的角度来看,衰退期的水体主要转化为裸地和林地,其比重分别占总土地面积的6.125%和3.721%. 水体退化为裸地的面积约为水体退化为林地的两倍,可见这一时期祁连山自然保护区出现了较为严重的水土流失现象. 因此,可以推断,水体面积的大幅度减少是这一时期生态安全指数下降的主要诱因之一.

(2)2010～2015年的涵养期. 图3与表6数据显示,2010～2015年间祁连山自然保护区土地生态安全指数变化极小,基本稳定. 说明这一时期保护区生态安全逐步开始修复,基本遏制住了水体的流失与裸地的增加,同时草地面积也有明显的增加,基本控制住了衰退期土地生态环境恶化的势头,效果显著.

(3)2015～2019年的修复期. 这一时间段保护区土地生态安全指数呈现出小幅波动中持续上升的态势,说明经过涵养期的生态修复,保护区在修复期实现了裸地的大幅度下降和水体的大面积增加,土地生态治理效果取得显著成效,土地生态环境表现出持续向好的发展趋势.

6 总 结

(1)本文以祁连山自然保护区土地利用类型的时空演变对土地生态安全状态的影响为研究重心,以土地利用转移矩阵和灰色关联分析方法,基于土地利用的净变化量、交换变化量和总变化量对祁连山自然保护区的土地利用变动情况分3个时期进行了深入的分析. 研究表明,土地利用类型的变动与祁连山自然保护区生态安全指数的波动存在较强关联关系,其中,水体面积的变动对土地生态安全的影响最为显著,同时,裸地面积和草地覆盖面积的变化也是影响这一地区土地生态安全的关键因素.

(2)本研究仅从土地利用类型变动的角度对祁连山地区的土地生态安全状态进行了探究,按照土地利用类型的大类与研究区的生态安全指数波动进行了宏观上的分析,可能存在着一定土地利用类型分类精度上的问题,因此,在进一步的研究中将关注土地利用类型面积及其时空变化的精准测度.

(3)研究结果显示,近年来祁连山自然保护区土地生态修复有了明显的成效,但与1985年相比仍然存在差距. 水源涵养依旧是祁连山自然保护区未来生态修复的主要工作方向. 应通过严格的水体治理,并加强林地用途管制,改善水文状况. 此外,还应该加大对石漠化沙化土地、工矿废弃地、生态重要区域的裸地的治理等措施,以提升保护区土地生态的自我修复能力.