王 南,仲 涛,郑 江 华,2,孟 乘 枫

(1.新疆大学 地理与遥感科学学院,新疆 乌鲁木齐 830017; 2.新疆大学 绿洲生态教育部重点实验室,新疆 乌鲁木齐 830046)

0 引 言

2021年8月,政府间气候变化专门委员会(IPCC)第六次评估(AR6)第一工作组报告指出,未来20 a地球表面变暖,温升将会达到1.5 ℃或1.6 ℃[1]。气温的升高加速了高海拔地区积雪的大量融化,补充了内陆干旱区的淡水资源,但可能引起冰川湖泊溃决洪水及泥石流灾害[2]。高分辨率遥感影像为冰湖监测及其潜在危害评价提供了支持。刘娟等[3]基于Landsat TM/OLI遥感影像资料,分析了帕隆藏布流域冰湖时空分布及其动态变化,并对潜在危险性冰湖进行了判别和评估;王翔等[4]基于野外调查和FLO-2D水文动力学模型等数据与方法,对藏西南地区典型冰湖强宗克措和吉莱普措开展动态变化分析、地貌调查和溃决泥石流模拟;King等[5]研究了20世纪70年代以来地表碎屑覆盖和冰川湖泊对喜马拉雅冰川质量损失的影响。

喜马拉雅山区是冰湖研究热点区域,而其他区域的冰湖研究则相对较少。当前,对昆仑山地区湖泊及冰川变化的研究主要有:沈永平等[6]在研究新疆冰川、积雪对气候变化的响应时发现昆仑山地区冰崩灾害随着气温升高而呈增加趋势;Bazai等[7]根据喀喇昆仑库尔多平冰川的涌动特征对该地区冰湖溃决形成的时间进行了估算;张太刚等[8]回顾亚洲高山区冰湖溃决洪水事件时,发现亚洲冰坝湖多分布于冰川资源丰富的天山山脉和喀喇昆仑山,未来几十年该地区冰湖溃决洪水频率将会增加。因此确定适用于昆仑山冰湖未来变化的预测模型,并讨论温度和降水等气候因子与预测结果之间的关系,将对该地区未来发展规划起到重要的作用。

元胞自动机(Cellular Automata,CA)是一种时间、空间、状态都离散,空间相互作用和时间因果关系皆局部的网格动力学模型[9],许多土地利用/覆盖预测模型是以此为基础而开发建立的。Wu[9]最先运用多标准评价模型,在CA模型中引入了多种约束因素,在用CA模型解决城市发展问题时考虑了宏观外部因素的影响;美国克拉克大学制图技术与地学分析实验室在此基础上,开发了微机地理信息和图像处理系统IDRISI(Idrisi for windows)[10]。国内学者进一步探索了CA模型的扩展以及构造规则。张显峰等[11]建立了城市土地利用/覆盖演化过程模拟预测LESP(Landuse Evolution Simulation land Prediction Model)模型,并对包头市的城市增长进行了模拟;Liu等[12]在神经网络模型基础上,结合自适应惯性竞争机制建立了未来土地利用/覆盖变化模拟模型FLUS(Future Land Use Simulation);Liang等[13]基于CA开发了新型土地利用/覆盖变化模拟模型PLUS(Patch-generating Land Use Simulation Model),能够有效处理人与自然共同因素影响下的土地利用/覆盖变化过程中的不确定性。

目前,土地利用/覆盖预测模型在模拟与预测城市扩展动态过程等方面应用较广,但是关于预测冰湖未来变化的研究较少,并且缺乏将已有预测模型应用到昆仑山冰湖未来变化的研究。因此,本文选取3种土地利用/覆盖预测模型(多准则CA-Markov模型、FLUS模型、PLUS模型),分析其在昆仑山冰湖预测中的适用性,从中选择出最适用于昆仑山冰湖预测的模型,预测昆仑山未来冰湖变化趋势,为该地区冰湖溃决、泥石流等灾害风险评估提供重要的参考依据。

1 研究区概况与数据处理

1.1 研究区概况

昆仑山是亚洲中部大山系,也是中国西部山系的主干。山脉西起帕米尔高原东部,全长约2 500 km,平均海拔超过4 000 m,宽130～200 km,总面积达50多万km2[14]。地势西高东低,地形以较为陡峭的山坡和较为高大的山峰为主,在北起伊尔克什坦山口直至新疆与西藏行政边界处存在较多山体与斜坡,东部地势相对平坦,冰湖在此处广泛分布(见图1)。夏季冰川和积雪大量消融是昆仑山地区冰湖水量补给的主要来源,因此该地区冰湖水量季节差异较大[15]。

图1 昆仑山地区冰湖地理位置Fig.1 Study area of glacier laks in Kunlun Mountain

1.2 数据来源及处理

本研究使用GEE(Google Earth Engine)平台获取昆仑山地区2000～2020年Landsat-5/7/8遥感影像数据,并采用随机森林法对影像进行监督分类。通过目视解译选取样本点,建立4种ROI(Regions of Interest):植被、裸地、冰雪、冰湖。样本点中70%用于分类,30%用于精度验证。通过构建混淆矩阵,计算得到总体分类精度(Overall Accuracy,OA)平均值为0.92,Kappa系数平均值为0.89。为减小由于季节差异带来的积雪融雪而造成的数据误差,使用云量较少、变化稳定的9～10月影像[16]。由于2000年、2004～2007年、2013年的影像云量过多,对数据质量带来的影响难以规避,所以在相应年份改为使用8～11月份的影像。

将监督分类结果与原始影像以及冰川、冰湖编目数据进行对比,通过人工目视解译进行冰川、冰湖提取。之后基于已有的关于冰湖提取方法[17],以冰川末端为中心建立10 km缓冲区,以实现对冰湖的二次精提取。

DEM数字高程数据采用数字高程数据产品SRTM DEMUTM,空间分辨率为90 m。该数据可用于提取研究区内的山体阴影、海拔等地形因子,进而对冰湖的分布进行空间分析。数据来源于地理空间数据云(http:∥www.gscloud.cn/)。

气象观测数据来自FLDAS数据集(https:∥disc.gsfc.nasa.gov/)1990～2020年的气温和降水数据,时间分辨率为1个月,空间分辨率为0.1°×0.1°。

冰川编目数据来自GLIMS全球冰川编目数据以及2017～2018年中国西北冰川编目数据[18],空间分辨率为30 m。冰湖编目数据来自高亚洲冰湖编目数据集(http:∥www.ncdc.ac.cn/)。

2 研究方法

2.1 多准则CA-Markov模型

基于IDRISI软件GIS Analysis下的CA-Markov模型将CA和Markov模型结合起来,运用转移概率矩阵来模拟土地覆盖格局随时间的变化[19]。首先利用IDRISI模型中的Markov模块,获得研究区2000～2010年和2010～2020年地物类型变化的转移概率矩阵:

(1)

土地覆盖类型转移适宜性图集指的是某一种地物类型转变为其他地物类型的概率图,适宜性图集是转换规则的一个重要组成部分[20]。本研究从高程、坡度、温度、降水4个方面考虑各因素对研究区土地覆盖变化情况的影响,应用IDRISI软件中MCE(多准则评价)模块的布尔相乘方法将各地物类型的约束条件进行合并,得到各种地物类型的适宜性图像,最后使用软件中的Collection Editor模块,按照各地物类型的顺序进行排列,最终合成研究区地物类型转变适宜性图集。

在确定P后,可用式(2)对未来土地覆盖类型进行预测。

S(t+1)=Pij×S(t)

(2)

式中:S(t),S(t+1)分别为t,t+1时刻的土地覆盖状态;Pij为土地覆盖类型转移概率矩阵。

在IDRISI软件CA-Markov模块下导入2000～2010年马尔科夫转移矩阵、2010年的昆仑山地区土地覆盖类型栅格图,以及制作的适宜性图集,循环次数设置为10,采用默认的摩尔型邻域5×5滤波器,获得2020年模拟影像。

2.2 FLUS模型

FLUS模型的原理是在利用人工神经网络模型算法(ANN)对基期土地覆盖数据和各驱动影响因子数据进行运算的基础上,估算出区域内各土地覆盖类型的发展概率[21]。之后在元胞自动机模块中,通过自适应惯性竞争机制可以定义自适应系数,该系数可以调整每类用地的惯性,并通过反复迭代未来不同情景下的土地覆盖量和当前各类用地的实际数量,最终模拟出未来特定情景下的土地覆盖分布情况[22]。自适应惯性系数定义式如下:

(3)

参考刘晓娟等[22]的研究方法,引入分块并行计算技术对昆仑山地区土地覆盖变化进行模拟,并行方法采用块状存储模式。采用人工神经网络整合并计算研究区域内每种土地覆盖类型在每个像元上的出现概率,获取研究区各类地物的适宜性概率,采样速率设置为20,隐藏层数为12。之后,将2010年土地覆盖类型图、驱动因素(高程、坡度、温度、降水)以及土地覆盖类型适宜性概率输入自适应惯性机制的元胞自动机模块,运行得到2020年土地覆盖类型模拟图。参考已有研究成果[21,23]并考虑昆仑山地区地物类型变化特征,依据研究区实际土地覆盖类型间的转换量对各个地类的邻域权重进行设置(裸地0.500、冰湖0.023、植被0.167、积雪0.310)。邻域权重参数范围为0～1,越接近1代表该用地类型的扩张能力越强[24]。

2.3 PLUS模型

(4)

式中:d的取值为0或1(若d=1,表示有其他土地利用类型转变为k类土地利用类型;当d=0,表示土地利用类型转变成了除k以外的其他土地利用类型);x是由若干驱动力因子组成的向量;函数I是决策树集的指示函数;hn(x)是向量x的第n个决策树的预测类型;M为决策树的总量。

最后结合CARS的随机斑块生成、过渡转移矩阵和阈值递减机制实现优化,确定最终土地覆盖类型分布情况[25]。

对2010～2020年昆仑山地区土地覆盖扩张数据和4类驱动因子(高程、坡度、温度、降水)进行LEAS用地扩张策略分析,生成每类用地发展概率。决策树数目20,采样率0.1,训练RF的特征个数为16,并行线程数量设置为4。之后基于2010年土地覆盖类型图与每类用地发展概率,使用CARS模拟得到2020年土地覆盖类型图。经过多次调试参数,根据各土地类型扩张面积的占比[13,26]将其扩张能力依次定义为:裸地0.500、冰湖0.023、植被0.167、积雪0.310。斑块生成阈值即递减阈值设置为0.9,扩散系数为0.1,随机种子最大比例为0.0001(取值0～1,数值越小越紧凑,反之则越分散)。

2.4 气候倾向率

气候倾向率可用来反映气候要素的变化趋势。一般是通过一元线性回归方程求得,即通过最小二乘法构建时间序列与气候要素之间的回归方程[19]。表达式为

yi=a+bxi(i=1,2,3,…,n)

(6)

(7)

式中:yi表示气候要素值;xj表示时间序列值;a为常数,一般用10b表示气候倾向率,单位为℃/10a或mm/10a。

3 结果与分析

3.1 模型精度检验

多准则CA-Markov模型、FLUS模型、PLUS模型的空间模拟情况如图2所示。PLUS模型在模拟各地物类型空间分布上优于多准则CA-Markov模型与FLUS模型,其模拟的2020年昆仑山地区土地覆盖类型空间分布与该地区实际情况一致性较好,但在模拟过程中存在裸地斑块丢失的情况;而FLUS模型则存在错误估算某类土地覆盖类型发展概率的问题,例如在昆仑山西南部本应有大量积雪覆盖,而FLUS模型模拟结果则显示该地区有部分植被分布。

图2 多准则CA-Markov模型、FLUS模型及PLUS模型空间模拟结果对比Fig.2 Comparison of space simulation results bymulti-criteria CA-Markov,FLUS and PLUS

为验证3种模型的模拟精度,将3种模型预测出的研究区2020年各类地物规模与2020年的实际值进行比较,计算Kappa 系数(见表1)。当Kappa≥0.75时,模拟结果和实际结果的一致性较好,精度较高;当0.40≤Kappa<0.75时,精度一般;当Kappa<0.40时,精度较差[27]。本研究采用多准则CA-Markov模型得到的预测结果总体精度(0.86)达到0.75以上,并且对研究区4种土地类型模拟的相对误差整体要小于另外两种模型;而采用FLUS、PLUS模型模拟的精度分别为0.63,0.67,均小于0.75。

表1 多准则CA-Markov模型、FLUS模型及PLUS模型模拟误差分析Tab.1 Error analysis of simulating on multi-criteria CA-Markov,FLUS and PLUS simulation

综合3种模型预测结果的Kappa系数与空间模拟情况,得出多准则CA-Markov模型预测结果与实际结果一致性较好,更符合昆仑山地区土地覆盖变化规律。以2015年昆仑山地区土地覆盖类型图为基础,以2010～2015年转移概率矩阵、制作的适宜性图集为转换规则,模拟2020年该地区土地覆盖变化,并对模拟结果进行精度检验。多准则CA-Markov模型以5 a、10 a为时间间隔模拟得到的Kappa系数均大于0.75(0.83、0.86),说明该模型模拟精度高,模拟预测的效果较好,可信度较高,适用于模拟研究区2025,2030年土地覆盖类型变化情况。

3.2 昆仑山冰湖变化特征分析

3.2.1冰湖变化时空特征分析

提取出的冰湖数据进行统计分析结果表明,2021年间昆仑山地区冰湖总体表现出数量增多、面积增大的趋势。2000年,该地区共有冰湖642个,冰湖面积共计28.03 km2;2020年,冰湖数量增长至894个,总面积达50.83 km2(见图3)。冰湖的总数量与总面积增加分别以平均每年2.22%和3.65%的速度增加。2015年之后,冰湖总数量与面积变化不明显,基本呈稳定增加趋势。

图3 2000～2020年昆仑山地区冰湖数量与面积变化统计Fig.3 Statistics on the numbers and areas of glacial lakes in Kunlun Mountain from 2000～2020

本研究在张连成等[15]的研究基础上,以83°E为界限将昆仑山地区分为东西两部分,分别讨论冰湖数量和面积的变化情况(见图4)。由图4可知,2020年,昆仑山西部冰湖数量占总量的73.60%,冰湖面积占全部面积的81.11%,研究区总体冰湖数量和面积表现出明显的西多东少的特点。西部冰湖的数量和面积增长剧烈,近20 a,西昆仑山冰湖平均以每年11.6个的速度增长;东部冰湖数量增长速度较慢,年均增长速度接近于0。东、西昆仑山冰湖面积变化速率较为接近,分别以3.78%、4.26%的年均增长率进行扩张。

图4 2000～2020年昆仑山地区东西部冰湖的面积与数量对比Fig.4 Comparison of the areas and numbers of glacial lakes in the east and west of the Kunlun Mountain in 2000～2020

3.2.2土地覆盖类型变化

利用IDRISI软件Markov过程获得的昆仑山地区2000～2010,2010～2015年和2015～2020年土地覆盖转移矩阵,分类汇总2000～2020年的研究区土地覆盖类型面积及变化情况(见图5)。

图5 2000～2020年土地覆盖转移桑基图Fig.5 Land cover transfer Sankey diagram from 2000 to 2020

昆仑山土地覆盖类型以植被、裸地、冰川和积雪、冰湖为主。由面积变化可知,昆仑山地区2000～2020年,冰川和积雪始终是该地区土地覆盖类型中面积增幅最大的地物类型,在2000～2010,2010～2015年和2015～2020年分别增加19 664.17,19 643.48 km2和12 867.44 km2;裸地面积呈持续减少趋势,减少量分别为19 548.59,21 119.68 km2和22 924.12 km2,主要转向植被、冰川和积雪;冰湖面积呈现持续增长趋势,在2000～2010,2010～2015年和2015～2020年按增长量分别为10.52,6.04 km2和5.27 km2,其中裸地、冰川和积雪向冰湖转化的面积较大,分别占冰湖转入总面积的42.39%～59.75%、25.69%～41.83%。

3.3 昆仑山冰湖未来变化趋势

通过对3种模型预测结果的Kappa系数与空间模拟情况的比较,得出多准则CA-Markov模型为最适用于昆仑山冰湖预测的模型。以2020年昆仑山地区实际土地覆盖数据、2015～2020年马尔科夫转移矩阵为基础,获得2025年该地区土地覆盖预测数据,同理获得2030年土地覆盖预测数据(见图6)。

经统计得到:2025年冰湖面积及数量预计达到54.16 km2和981个,与2020年相比分别增加了6.54%、9.73%;2030年冰湖面积及数量预计达到58.25 km2和1 042个,与2020年相比分别增加了14.59%、16.55%,与2025年相比增加了7.55%、6.22%。

将研究区以83°E为界分为东西两部分,进一步分析得到:2025年西昆仑山冰湖面积预计为44.22 km2,东昆仑山冰湖面积预计为9.94 km2;2025年东部冰湖面积与2020年相比增加约3.53%,而西部冰湖面积则将以7.25%的速率增长,由此可见,虽然研究区的大面积冰湖集中分布于东部,但2020～2025年面积增长速度较快的冰湖仍主要来自西部;到2030年,东昆仑山冰湖面积及数量增长为10.64 km2和272个,而西昆仑山冰湖面积及数量将达到47.60 km2、770个,与2025年冰湖面积相比,研究区西部、东部冰湖面积分别增加7.66%与7.09%,数量分别增加6.65%、5.02%,东部冰湖增长速度有所提升,且东西部冰湖扩张的速度差距有所减小(见表2)。

表2 2020、2025年及2030年昆仑山地区东、西部冰湖面积与数量Tab.2 The area and number of glacial lake in the east and west of the Kunlun Mountain in 2020,2025 and 2030

气温、降水关系到冰川生长退缩、地表径流大小、蒸发速度快慢,其变化还会影响冰湖的补给水源、湖盆水量以及湖盆结构[28]。因此,分析昆仑山地区气温、降水变化情况对研究冰湖的演化发展有重要意义。通过计算气候倾向率得到2000～2020年昆仑山西部降水倾向率为0.29 mm/10 a,气温倾向率为0.04 ℃/10 a;昆仑山东部降水倾向率为-4.92 mm/10 a,气温倾向率为0.16 ℃/10 a。研究区西部降水量增幅较小且温度整体保持平稳,该区域有大量冰川和积雪覆盖,在水热条件及积雪消融共同作用下,冰湖面积与数量将继续增加;研究区东部降水呈下降趋势,但由于气温增幅较西部地区更为明显,由此带来的冰川融水补充了一部分由于降水量减少而缺失的冰湖补给,因此东西部冰湖扩张的速度差距将有所减小。

在地形方面,昆仑山东部地势较低,坡度较为和缓,有利于大型湖泊的发育;西部地势较高且地形复杂,冰湖在此常呈小规模分布;而在昆仑山脉的北坡,是最干旱的亚洲大陆中心,山前年降水量小于100 mm,随山地海拔增高而略增,气候干燥,再加上昆仑山中段主脉向南略呈弧形,预计2025年、2030年内,由于干燥的大陆气候的影响,北部的狭长地带仍将被裸地所占据,冰湖以及积雪等其他地物类型在此则鲜有分布。

3.4 分析与讨论

综合3种模型预测结果的Kappa系数与空间模拟情况可以看出,多准则CA-Markov模型模拟结果的面积较为符合2020年真实的分类结果。陈柯欣等[29]研究发现多准则CA-Markov模型Proportional error参数通过影响转移概率矩阵,进而对模拟数量精度产生影响,该参数为0.15时,影像的分类精度一般可达到85%以上。本文在该模型Markov过程中将Proportional error参数设置为0.15,结果显示数量模拟效果较好;蒋小芳等[30]提出PLUS模型模拟空间位置分布的效果优于其他模型,且在越精细的空间尺度下的模拟精度越高。本研究采用90m空间尺度,该尺度下的各土地覆盖类型和驱动因子相关性不高,可能是模拟过程中部分裸地斑块丢失的原因之一[31];FLUS模型为了更好地模拟土地覆盖的动态性和不确定性,具有较低概率的地物类型也有机会被分配到目标地物类型[21],因此有可能出现本文中FLUS模型估算昆仑山地区土地覆盖类型发展概率时出现偏差的情况。

已有研究表明,气温上升是昆仑山地区冰川退缩的主要原因[32]。东昆仑山冰川规模较小,但冰川退缩速度较快,气温上升使得冰川消融大于降水增加引起的冰川积累,冰湖面积及数量因此保持增长趋势。青藏高原冰川响应气温变化的滞后时间为10 a左右[33],因此未来10 a内,该地区冰川将出现持续性的退缩。本文通过多准则CA-Markov模型获得2025、2030年土地覆盖预测数据,并通过进一步分析得到研究区东部降水呈下降趋势,但由于气温增幅较西部地区更为明显,由此带来的冰川融水补充了一部分由于降水量减少而缺失的冰湖补给,因此东西部冰湖扩张的速度差距将有所减小,该结论与前人研究昆仑山冰川、冰湖变化的结论一致。

本文在使用3种预测模型对昆仑山地区土地覆盖情况进行模拟过程中仍存在不足之处:仅基于经验与土地覆盖变化的历史特征进行参数设置,没有考虑未来不同情境下气候因素对土地覆盖类型造成的影响,因此如何构建更为全面、科学的预测模型仍需要进一步探讨。另外,依靠人工目视解译得到的监督分类结果即使经过了与其他数据集比对后的二次提取,依然存在与实际冰湖分布情况之间的误差。从水体和积雪指数入手探索精度更高的冰湖提取方法将是未来工作的重点。

4 结 论

本文基于GEE平台进行监督分类,分析了昆仑山地区冰湖2000～2020年变化时空特征。通过比较3种预测模型模拟结果,确定出适用于昆仑山冰湖的预测模型,并对2025、2030年冰湖可能的演变发展情况做出预测,得出如下结论:

(1) 20 a间昆仑山地区冰湖总体表现出数量增多、面积增大的趋势。2000年,该地区共有冰湖642个,冰湖面积共计28.03 km2;2020年,冰湖数量增长至894个,总面积达50.83 km2。研究区总体冰湖数量和面积表现出明显的西多东少的特点。2020年,昆仑山西部冰湖数量占总量的73.60%,冰湖面积占全部的81.11%,该区域冰湖的数量和面积增长剧烈,是研究区冰湖数量和面积增加的主要成分。

(2) 通过比较多准则CA-Markov、FLUS、PLUS模型预测结果,确定多准则CA-Markov模型为最适用于昆仑山冰湖的预测模型。该模型以10 a为时间间隔模拟得到的Kappa系数为0.86,并且对研究区4种地物类型模拟的相对误差整体要小于另外两种模型;采用FLUS、PLUS模型模拟的精度分别为63.46%、67.47%,均小于75.00%。

(3) 2025年冰湖面积及数量预计达到54.16 km2、981个,2030年冰湖面积及数量预计达到58.24 km2、1 042个。2025年西昆仑山冰湖面积预计为44.22 km2,东昆仑山冰湖面积预计为9.94 km2,东、西部冰湖面积增长速度分别为3.53%、7.25%,由此可见2020～2025年面积增长速度较快的冰湖仍主要来自西部;到2030年,东昆仑山冰湖面积及数量增长为10.64 km2、272个,而西昆仑山冰湖面积及数量将达到47.60 km2、770个,与2025年冰湖面积相比,研究区西部、东部冰湖面积分别增加7.66%与7.09%,数量分别增加6.65%、5.02%,东部冰湖增长速度有所提升,且东西部冰湖扩张的速度差距将有所减小。