珠峰北麓朋曲流域土地覆被变化及未来情景模拟

李存秀刘峰贵周源涛陈 琼周 强何文鑫田万霞

(1.青海师范大学 地理科学学院,青海 西宁810008;2.高原科学与可持续发展研究院,青海 西宁810008)

土地覆被/土地利用作为人类活动响应气候变化的重要途径之一[1],在全球变化中扮演着重要的角色,是全球变化的重要驱动力。工业革命以来,随着经济和科学技术的迅速发展,以及在全球变化的高速推动下,人地矛盾日益尖锐[2]。近年来中国特别强调,保障生态文明有效建设和可持续发展的关键是维持人地关系协调发展[3]。在生态环境脆弱、自然环境复杂的西藏地区,人地关系协调演进无疑是一项重大挑战。随着人地关系研究的不断推进,西藏地区人地关系优质协调发展的势头良好[4]。珠峰北麓朋曲流域是中国珠穆朗玛峰自然保护区的核心腹地,高山草甸作为环珠峰地区的主要植被类型,20世纪末至21世纪初,随着牧业比重的增大,环珠峰地区过牧现象严重,草地退化加剧,森林资源也处于明显的退化状态,退化面积大于恢复面积,区域生态环境遭到破坏,人地关系冲突日益明显[5],因此,加强环珠峰地区生态环境保护问题被提到重要地位。

从20世纪60年代第一次展开珠峰科学考察开始,珠峰自然保护区就成为高原人地关系的研究热点,研究聚焦人与自然相互作用的评估[6]、气候变化对环珠峰地区植被动态变化规律及影响等研究[7-8]。研究结果表明,自2002年以来,随着西藏退耕还林、还草工程的实施以及环珠峰地区生态保护工程的实施,珠峰南坡植被类型恢复良好,珠峰北坡植被稳定发展[9],但耕地面积明显减少。环珠峰地区土地覆被研究主要集中于植被现状和恢复研究,但对未来情景下的土地覆被动态的研究不足。

当前用于土地覆被变化模拟预测的模型主要有:CLUE-S模型、Logistic回归模型、Markov模型和FLUS模型等[10]。FLUS模型是对传统预测模型的改进,能有效解决传统预测模型中出现的参数设置和类型之间转换困难等问题,有效处理人为和自然因素双重影响下的土地覆被/土地利用变化不确定性;Markov模型对预测未来情景下的土地覆被/土地利用像元数量具有较高的实用性,但难以对变化进行空间表达[11]。为此,基于Markov-FLUS模型中和了Markov模型的缺点,以达到像元数量预测和空间表达的双重优势,适用于中小流域未来土地覆被情况变化的研究。

在青藏高原生态保护和可持续发展的大背景下,对耕地和生态安全保护将对未来情境下的流域水土响应与人地关系的和谐发展产生重要影响。因此,针对流域内土地覆被变化未来发展的研究非常之必要,既可以对保障粮食安全和保护生态安全提供建设性意见,又为下一步生态文明发展提供基础理论支撑。朋曲流域作为环珠峰北麓地区的人口中心,是珠穆朗玛峰国家级自然保护区的核心区与人地关系相对复杂的地带,同时作为珠峰地区重要的水文基本单元,统筹流域的相关研究,是构建国家生态安全屏障的重要研究对象。对朋曲流域未来土地覆被情况的研究,可以反映高海拔高脆弱性的特殊自然环境下未来土地利用情况的动态变化规律,以及不同土地利用限制情景下,土地覆被在未来近期和中期可能发生的演变。本文基于Markov-FLUS模型,运用不同情景约束因子,对朋曲流域短期(2030年)和中远期(2050年)的土地覆被类型进行模拟预测分析,对比两种类型下的耕地发展和生态保护发展,为促进朋曲流域耕地保护工作和实现生态空间的科学布局提供思路,为实现高原地区和谐的人地关系服务。

1 研究区概况

朋曲发源于喜马拉雅中段的希夏邦马峰,源头以冰川融水为主,大气降水补给为辅,河流自西向东在定结县境内折向南,切穿喜马拉雅山脉入尼泊尔,最后注入印度恒河,上游中国境内流域范围涉及西藏日喀则市南部的岗巴县、定日县、定结县、聂拉木县和萨迦县,全长376 km,流域总面积2.50×104km2,平均海拔5 259 m,支流主要有扎嘎曲、叶如藏布江、西宁藏布和卡得藏布等[12]。喜马拉雅中段的珠峰地区是青藏高原人类活动历史悠久且较为频繁的区域,该河流属于跨境河流,是喜马拉雅中段人类南北迁徙和贸易的主要通道,也是珠穆朗玛峰国家自然保护区的核心区域,因此,该流域受到广泛关注。

朋曲流域焚风效应明显,印度洋季风暖湿气流在喜马拉雅山脉南坡形成大量降水,越过山脊到达北坡后下沉,导致该地区常年气温较低,降水偏少,年均降雨量仅仅维持在600 mm左右,属于典型的高原大陆性气候,寒冷而干旱,由此北坡植被多为草甸和灌丛。位于珠峰南坡的少部分区域,受印度洋季风影响,降雨充沛[13],气候温暖湿润,形成大片原始森林。朋曲水系包括热曲藏布、扎嘎曲、罗别藏布、彭作浦曲、西纠曲、鲁鲁曲等小河流和主河道朋曲河[14]。据统计,2020年朋曲流域总人口为16.5万人,人口密度约为4.7人/km2,半农半牧,耕地主要分布在河谷地带,农作为以种植青稞、小麦、油菜为主[15]。

2 数据来源与方法

2.1 数据来源与处理

本文数据包括朋曲流域县界数据,空间驱动因子数据集和土地覆被数据集,研究区县界数据来源于2019年西藏自治区县级行政边界提取的数据;空间驱动因子数据集包含自然因子、社会经济因子和距离因子等数据;2010,2015和2020年3期土地覆被数据来源“地球大数据工程数据共享服务系统”[16](表1)。本文以2020年最新研发的高精度、全球性土地分类数据(GLC_FCS30-2020)为基准数据,该数据集共包含29个地表覆盖类型,更新周期为5 a。本文通过筛选所需要的土地覆被类型,将29个地表类型进行重分类操作,最后提取需要的5个土地覆被类型。数据来源详见表1。

表1 研究所用各类数据来源

2.2 FLUS模型

FLUS模型是主要用于当前和未来情景下的土地覆被模拟。该模型的原理源自元胞自动(selfadaptive inertia and competition mechanism CA),并且在其基础上做了进一步修改,增加了人工神经网络模型算法(ANN-based probability-of-occurrence estimation)和马尔科夫链计算(Markov chain)两大模块[18],模型的改善和发展进一步耦合了人类活动和自然要素之间相互竞争和转化的模拟优势[19]。

2.2.1 适宜性概率统计 人工神经网络模型可用以表征未来情景土地覆被类型和多项空间驱动因子之间的对应关系,将初始年份的土地覆被数据和各项驱动因子作为输入数据源,采用神经网络算法[20](ANN)整合并计算朋曲流域内5种土地覆被类型出现在每个像元上的概率[21],计算公式为:

式中:Nj(p,q)为隐藏层神经元j接收的输入层信号;wij为训练中输入层和隐藏层之间设定的权重[22]。

2.2.2 马尔科夫链预测未来像元数 利用马尔科夫链模块预测未来情景下的土地覆被类型像元数量[23],在2010,2015和2020年的数据基础上,预测2020年和2025年的土地覆被类型像元数量(表2),计算公式为:

表2 研究区2010—2050年多期土地覆被类型像元总数

式中:Tn+1为n+1时刻土地覆被状态;Tn为n时刻土地覆被状态;pij为各地类的转变适宜性概率。

2.2.3 领域因子权重计算 领域因子用以表征不同土地覆被类型的扩张能力[18],扩张能力范围为[0,1],范围数值大小与扩张强度成正比,当扩张数值越大接近于1时,表示该土地覆被类型在其领域范围内的扩张能力越强。当前用于确定权重参数的方法主要有:因子分析法、Moore邻域模型法[24]、无量纲化。无量纲化方法将原始数据进行无量纲化处理,根据每个指标值的差异,客观确定指标权重[25]。本文考虑土地类型扩张规律的客观性,使用朋曲流域各土地类型历史总面积变化数量进行无量纲化处理,来确定各地类的领域因子权重值[26](表3)。

表3 朋曲流域无量纲化得到的领域因子权重值

2.2.4 转换成本计算 转换成本反映各土地类型之间的变化规则,通常用矩阵表示,0表示两种地类之间转换成本高,不容易转换,1表示两地类之间转换成本最低,最容易转换。当前关于计算转换成本矩阵的方法主要有专家意见咨询、经验转换和政策导向转换规则,本文利用朋曲流域2010—2020年土地类型转移矩阵为基础,以专家咨询意见为辅,确定各地类的转换成本系数。

2.3 模型精度验证

本文将验证模型模拟精度高低的评价指标分为两类:分类别精度评价指标〔生产者精度PA(producer’s accuracy)〕和用户精度UA(user’s accuracy)检验[27])和综合精度评价指标〔总体精度OA(overall accuracy)和kappa系数(kappa coefficient)检验〕。其中,总体精度OA和Kappa系数主要检验模拟整体结果与实际结果之间的一致性,用以评价模拟图像的精度问题,两个结果的检测值越接近于1,表明两个图像之间具有较高的相似性[28],本文主要采用该方法进行模拟精度的验证。以朋曲流域2015年为初始年份数据,对影响因子归一化处理,采用FLUS模型的马尔科夫链模块预测2020年的土地类型数量,计算得到朋曲地区2015—2020年的土地覆被转移适宜性概率数据,然后再对2020年的朋曲流域土地覆被进行模拟,将模拟结果与2020年实测的土地覆被类型进行精度验证,利用FLUS模型模拟出的2020年数据结果总体精度(OA)为90%,kappa系数为0.82。其中草地的分类精度最优,PA和UA均在90%以上;其他用地类型的分类精度次之,PA和UA均超过88%;林地的UA高达93%,但是PA仅有66%;水域的分类精度相较于前3种类型较低,PA和UA仅有65%和77%;耕地的分类精度是最低的,PA和UA均不超过60%,分类精度一般(图1)。总体而言,FLUS模型对朋曲流域土地覆被类型的模拟结果精度较高,达到预期目标,模型可以用于本文关于未来情景土地覆被模拟预测研究。

图1 朋曲流域2020年土地覆被现状与模拟结果比较

3 结果与分析

3.1 历史时期朋曲流域土地覆被变化

根据2010—2020年朋曲流域土地覆被面积变化表和空间分布图(图2)可知,从2010—2020年期间,林地的面积变化最高,其次是草地,再次是耕地,其他用地和水域的面积变化是最小的。林地主要分布在聂拉木县、定结县和定日县南部地区,其面积变化最为明显,10 a间增加了778.16 km2,主要增加的区域为聂拉木西北部和定日县中部地区;其次是草地,草地是研究区内面积最大的地物类型,广泛分布于研究区内,面积呈显著较大的减少变化趋势,10 a以来共减少726.89 km2,草地面积减少区域主要为林地的扩张;再次为耕地面积变化,河谷地区零星呈带状分布的耕地,面积呈减少趋势,但是减少面积较小,为67.78 km2;其他用地面积稍有增加,为8.46 km2;朋曲流域的河流受到冰川融雪和大气降水的多重补给,常年径流量稳定,水域面积基本无变化;造成上述这种变化的主要原因是朋曲流域是以牧业为主农业为辅的地区,畜牧业是当地农牧民最重要的经济来源,草地面积在该流域占很重要的份额(表4)。随着气候的变化和社会经济的发展,越来越多的农牧民开始脱离草场,从事非农牧活动[14],人类活动区域逐渐向城镇集中,耕地面积逐渐减少,同时随国家“退耕还林”与珠峰国家级自然保护区的建立,撂荒地逐渐被林地取代。

图2 朋曲流域2010—2020年土地覆被空间分布

表4 朋曲流域2010—2020年土地覆被面积变化

3.2 未来情景土地覆被变化分析

利用FLUS模型模拟预测未来短期(2030年)和未来中远期(2050年)朋曲流域土地覆被状况,但由于马尔科夫模块预测的土地覆被类型像元数量结果,可分为3种限制情景的模拟预测,未来情景的土地覆被类型预测需要同时考虑不同限制因子对土地扩张的约束。因此,设置多个约束用地变化的限制因子(表1),对未来短期与未来中远期不同限制情景下的流域土地覆被变化进行模拟预测,其结果如下:

3.2.1 自由发展情景 自由发展情景旨在不考虑相关约束用地发展的政策因子,仅考虑在自然和人为影响下土地覆被的变化情况,图3和图4反映了在此情境下未来短期(2030年)和未来中远期(2050年)朋曲流域土地覆被空间分布状况与面积变化,较2020年:①未来短期土地覆被变化中,林地和其他用地面积呈增加趋势,而耕地、草地和水域较历史时期呈减少趋势,其中,其他用地的面积变化是所有土地覆被类型中变化最为明显的类型;②未来中远期土地覆被变化相较于2020年有些地类有明显的向外扩张趋势,林地面积的变化最为明显,较2020年增加3 498.35 km2,从空间分布图上看出2050年林地面积扩张明显,扩张范围主要集中在岗巴县和定结县范围内,同时陈颖峰等[29]人的研究结果也表示,在2000—2009年间,珠峰东南区域(岗巴县和定结县)林地等植

图3 朋曲流域自由发展情景下土地覆被类型面积变化比较

图4 朋曲流域自由发展情景下土地覆被类型面积变化

被类型增加明显,研究结果与以上结论吻合。在不受政策等限制因子去约束土地类型扩张的条件下,其他用地类型将快速扩张,扩张的面积主要占据了草地、耕地和水域等用地类型,若不加以相关政策因子限制,将直接威胁研究区内的生态环境和粮食安全。

3.2.2 耕地保护情景 基于耕地保护情景的模拟预测,以基本农田保护区为限制性因素,在该情景下的耕地发展稳定性较强,耕地类型向其他类型转化的条件受到限制。未来短期和未来中远期的土地覆被变化在此限制情景下发生了较大的变化(图5—6)。①在此情景下,未来近期耕地面积126.8 km2,较2020年耕地面积无明显变化,耕地减少的速率变缓,耕地保护的措施得到有效发挥。面积变化最明显的是其他用地类型,较2020年面积增加393.19 km2,这部分增加的面积主要来自于对草地面积的转换,致使草地面积减少幅度较大。朋曲流域的土地覆被类型复杂,其他用地类型里面包含的二级类型又充满了复杂性,其他用地类型的扩张,一方面有朋曲流域少部分的建设用地和冰川积雪地向外扩张的原因,而最重要的一方面主要为草地转化为裸土地,成为其他用地外扩的主要原因;②未来中远期林地的面积变化最大,较2020年面积增加3 485.55 km2,平均每年增加116.2 km2,这部分增加的面积主要来源于对草地和其他用地类型的转化。耕地变化是其余地类中面积变化最小的一类,耕地保护情景约束了未来中远期耕地向其他地类转出的趋势。从面积变化趋势分析,加入农田保护区作为限制因子的朋曲流域土地覆被变化,减缓了耕地转化的速率,可以有效保护耕地。

图5 朋曲流域耕地保护情景下土地覆被类型面积变化比较

图6 朋曲流域耕地保护情景下土地覆被类型面积变化

3.2.3 生态保护情景 朋曲流域大部分区域处于珠穆朗玛峰自然保护区范围内,生态保护是该区域可持续发展的基本方向。在生态保护情景下,从各用地类型数量变化见图7—8。①此情景下未来短期土地覆被类型变化较大的仍集中在草地和其他用地类型,相较于2020年的草地发展,未来近期草地面积减少362.67 km2,其他用地面积增加389.29 km2。此情景下的土地覆被类型面积变化与耕地保护情景下的面积变化类似,但是在空间分布上却出现了较大的不同。受自然保护区的影响,林地范围在研究区东南方向扩张迅速,并且呈大块连片分布;②未来中远期覆被类型面积变化也是主要草地、林地和其他用地,林地仍然是变化最大最为明显的地类,而草地和其他用地仍然呈减少趋势,面积受到进一步压缩,草地和其他用地的转出的主要方向为林地。但值得一提的是,此时期水域面积较2020年增加了117.538 km2,由于对气候变化的响应,珠峰自然保护区变暖趋势要比中国其他区域要更加显著[30],致使珠峰地区湖泊冰雪融水补给量增加,面积扩大,这与王毅等[31]人的研究结果相吻合。

图7 朋曲流域生态保护情景下土地覆被类型面积变化比较

图8 朋曲流域生态保护情景下土地覆被类型面积变化

生态保护情景下,土地覆被变化特征主要为草地的减少以及其他用地类型的增加,耕地面积变化在整体变化中所占比例不高,但依旧呈减少趋势。此情景下土地覆被类型变化较大主要集中在草地和其他用地类型,但是相较于自由发展情景和耕地保护情景下的土地扩张,生态保护下的土地覆被总体而言变化最小,土地转化速率趋于放缓,这对保护生物多样性,保障生态安全区域生态安全具有重要价值,同时我们也关注到,其他利用地与草地在此情景下变化仍最明显,究其原因可能与模型中使用的马尔可夫链有关。执行农田保护区和自然保护区的同时限制的相关政策,将对本区域的生态恢复、发展和经济建设提供重要的参考,有助于实现区域的可持续发展。

4 讨论与结论

4.1 结论

(1)利用Markov-FLUS模型,模拟2020年朋曲流域的土地覆被状况,和2020年实际的土地覆被现状进行对比。经过验证,模型分类的总体精度(OA)为90%,kappa系数为0.82,精度较高,达到预期模拟精度目标,所以可以采用此模型进行未来情景的土地覆被模拟预测。

(2)通过统计分析2010,2015和2020年土地覆被状况,林地和草地的面积变化最为明显,草地面积在10 a间减少726.89 km2,林地面积增加778.16 km2,主要是由草地转换为林地,主要是因为西藏地区近些年的植树造林政策的有效实施。

(3)在3种情境下,未来近期(2030年)和未来中期(2050年)的土地覆被均发生明显变化。未来近期和未来中期的土地覆被变化分别在连续3种情境下的表现一致,未来短期土地覆被变化最明显均集中于草地和其他用地两个类型,未来中远期覆被面积变化最明显地集中在草地、林地和其他用地。也就是说,在随着气候变化、经济发展和政策的逐步推行之下,朋曲流域的林地面积将迅速扩大,而这部分的转入是草地和其他用地的转出获得。

4.2 讨论

不同发展情景下的土地覆被变化研究对区域的经济和可持续发展具有特别重要的影响,也对当前朋曲流域的土地覆被空间格局变化具有良好的指导意义。朋曲流域作为珠峰自然保护区的主体部分,既要保障流域内粮食生产安全,又要保证流域内生物多样性的有序发展。尽管流域内的人口比重不大,当前的粮食生产能满足当前的人口规模,但是未来人口的不断增长会在一定程度上加重对粮食需求的压力,所以耕地保护情景下的土地覆被变化研究具有非常之必要性。朋曲流域生态保护作为珠峰自然保护区的核心问题,生态意义重大,生态保护情景下的朋曲流域土地覆被研究可以进一步弥补生态建设基于现状分析上的短板问题。虽然关于朋曲流域未来短期(2030年)和未来中远期(2050年)的模拟研究成果较好地拟合了未来土地发展状况,但是仍存在一些不足:①在使用模型中,未利用地与草地始终是变化最明显的两种覆被类型,这可能与模型中的马尔可夫链的特性有关,马尔可夫链的引用虽解决了随时间变化的土地覆被类型的预估问题,但作为一种基于概率论和数理统计的模型,在受全球变化与人类社会发展等一系列复杂机制影响的土地利用变化研究中,模型的预测效力还需进一步讨论;②模型验证中,实测结果与模拟结果kappa值为82%,虽满足了模型要求的75%以上的精度,但对于小流域的预测来说,仍可能与实际产生较大的差距。这种偏差的产生可能与青藏高原特殊的自然环境条件和区域政策有关,模型针对青藏高原进行区域性的改进和更多驱动因子选择与其权重的研究将是下一步的工作方向;③FLUS模型在进行未来情景模拟预测时,所选取的驱动因子指标具有一定的主观性,驱动因子指标的选取多参照前人研究成果中所选取的指标,因为所选取的研究区不同,受地势、海拔、气温、降水、人口、经济、宗教等因素的影响,不同研究区内对于土地利用的方式不同,不同地类发展的趋向也就各不相同,这就会在之后的模拟研究中增加模拟误差,模拟结果与实际发展方向背道而驰。另外在选择驱动因子时,无法准确分析驱动因子选取的多少是否对模拟结果产生重要影响,一系列问题还需进一步模拟验证。