1、田梦祺,段克勤,石培宏.基于 Google Earth Engine 平台的青藏高原冰川变化研究以普若岗日冰原为例 J.地理科学,2023,43(6):943-951.Tian Mengqi,DuanKeqin,Shi Peihong.Glacier changes on the Qinghai-Tibet Plateau based on Google Earth Engine:A case study of the Purog Kangri glacier.Scientia Geographica Sin-ica,2023,43(6):943-951.doi:10.13249/ki.sgs
2、.2023.06.001基于 Google Earth Engine 平台的青藏高原冰川变化研究以普若岗日冰原为例田梦祺，段克勤，石培宏（陕西师范大学地理科学与旅游学院，陕西 西安 710119）摘要摘要：基于 Google Earth Engine（GEE）云计算平台，以青藏高原中部的普若岗日冰原为例，探讨了快速批量处理遥感数据并获取冰川面积变化的方法。结果表明，利用 GEE 平台结合像元级影像合成和最小归一化积雪指数（Nor-malized Difference Snow Index,NDSI)，可快速获取年度合成影像，提取冰川边界。通过对 19882020 年共 512 景遥感影像合成计
3、算，获得了 19882020 年普若岗日冰川面积变化特征，分析发现：33 a 间冰川面积减小了 31.08km2，年均退缩率为 0.23%/a，冰川面积退缩主要发生在冰川末端，部分末端在 20 a 间退缩了 1 506 m。相对于冰川的实地观测或传统遥感监测法，本方法提高了影像利用率和数据处理效率，可快速实现冰川面积变化长时序分析，为现代冰川变化研究提供了新途径。关键词关键词：Google Earth Engine；普若岗日冰原；云计算中图分类号中图分类号：P343.6文献标识码文献标识码：A文章编号文章编号：1000-0690(2023)06-0943-09 在全球变暖背景下，青藏高原冰川正
4、在发生快速的变化1-3，导致冰川处于失稳状态中4，冰川的退缩已经引起了局地水循环的改变5，并引发一系列与冰川退缩相关的灾害6。因此，快速明确青藏高原冰川的变化趋势及幅度，有助于丰富和理解气候变化背景下青藏高原冰川失衡失稳的科学内涵、过程机理，对减缓和应对这一演化趋势具有十分重要的科学意义和应用价值。青藏高原冰川因其高海拔和难以到达的特点，对冰川长期且大规模的现场观测受到极大限制，学者只在容易到达的少数冰川进行了长期的监测7-8。中国第一、二次冰川编目项目，通过野外实地考察或遥感影像目视解译获取冰川变化信息，工作耗时长、效率低9。随着卫星遥感技术的发展，遥感监测逐渐成为冰川研究的重要方法10。利
5、用各类冰川提取算法结合目视解译，可获得青藏高原不同山系冰川分布、冰川面积及长度的变化11-12。通过对比不同时期卫星资料，发现青藏高原大多数冰川呈现加速消融的状态13-15。然而，青藏高原冰川变化存在显著的时空异质性，为明晰冰川变化的气候响应差异机理，首先需要通过海量数据获取大范围冰川变化时间序列数据。但目视解译遥感图像存在工作量大、周期长和主观判断误差等问题，基于单期影像或人工拼接的计算机自动解译，难以快速实现海量数据的同步计算16，在大尺度长序列的研究中耗时较长，导致大多数研究只是集中在个别山系和某一特定时间段。因此，为快速获取青藏高原冰川变化数据，亟需发展一种快速批处理方法。目前冰川遥感
6、监测算法主要包括基于波段组 收稿日期收稿日期：2022-01-22；修订日期修订日期：2022-04-07基金项目基金项目：第二次青藏高原综合科学考察研究项目(2019QZKK0201)、国家自然科学基金项目(41571062)、中国博士后基金项目(2017M610622)资助。Foundation:The Second Tibetan Plateau Scientific Expedition and Research(STEP)Program(2019QZKK0201),Nation-al Natural Science Foundation of China(41571062),the
7、China Postdoctoral Science Foundation(2017M610622).作者简介作者简介：田梦祺（1998），女，山东潍坊人，硕士研究生，主要从事冰川变化研究。E-mail:通信作者通信作者：段克勤。E-mail: 第 43 卷第 6 期地理科学Vol.43 No.62023 年 06 月Scientia Geographica SinicaJune,2023合的分类方法、传统机器学习分类方法和基于神经网络的深度学习方法。波段组合分类算法通过对冰川的特征波段进行组合及运算，自动或人工确定阈值提取冰川边界，如利用归一化积雪指数（Normal-ized Differe
8、nce Snow Index,NDSI)和波段比值阈值法等。波段组合法在提取冰川范围时易受到积雪、地形及阴影等因素影响17，通常需要结合目视解译提高冰川识别的准确率18。在冰川表碛覆盖区域，由于其光谱特征与非冰川区类似19，并且缺乏连续、高质量的影像，因此表碛覆盖型冰川的提取一直是冰川分类研究中的难题。NDSI 积雪指数等传统分类方法在洁净冰川区识别效果较好，但不能自动准确地提取表碛覆盖区的冰川20，对此，利用随机森林的机器学习算法可较高精度的识别喀喇昆仑山区的冰上碎屑（即表碛覆盖）区域21。遥感大数据时代的到来，数据量的爆发式增长，促使新的遥感数据处理平台，即云计算平台的产生和发展。早期冰川
9、变化遥感监测主要基于本地平台，需要下载大量影像，耗费人力和网络资源。遥感云计算平台的出现，解决了数据获取困难、本地系统建设成本高等问题，使快速获取并处理空间大范围冰川变化数据成为可能。Google Earth Engine(GEE)是目前应用最为广泛的遥感云计算平台，其强大的数据存储及处理分析能力，改变了传统的单机遥感影像处理模式，拓宽了遥感数据的应用范围和时空尺度22，极大提高了遥感监测研究的数据处理效率。GEE 平台在土地利用变化23、森林植被动态监测24、水体时空变化25-26等领域已取得广泛应用。目前也有个别研究基于 GEE 平台对高亚洲冰川变化进行了分析27，但提取方法仍在探索之中。
10、本文的目的是基于 GEE 云计算平台，以青藏高原腹地的普若岗日冰原为例，结合最小 NDSI 合成算法提取冰川边界，探讨传统波段组合算法与云平台结合进行长时序冰川变化动态监测的可行性，旨在为快速获取冰川变化特征提供一种新途径。1 数据与方法普若岗日冰原（8920E8950E，3344N3404 N）位于青藏高原中部（图 1），是中低纬度面积最大的冰原。普若岗日冰原由数个冰帽型冰川组成，20 世纪 80 年代冰川总面积约为 422 km2，冰川雪线海拔 5 6205 860 m，最高峰海拔 6 482 m 28。1.1 数据数据研究采用的影像数据为 19882020 年（由于Landsat 7 E
11、TM+影像存在严重的条带问题，故未使用 2012 年的数据）经过辐射校正和大气校正的Landsat 5 TM 和 Landsat 8 OLI 地表反射率数据集，空间分辨率为 30 m，时间分辨率为 16 d。为确保影像能够覆盖整个研究区，使用了条带号为 139 和140，行编号为 36 和 37，共计 512 景影像（https:/ JavaScript API 接口在线访问和云端处理。图 1 普若岗日冰原地理位置Fig.1 Location of the Purog Kangri glacier944地理科学43 卷 选用 6.0 版本的 RGI（the Randolph Glacier I
12、n-ventory）冰川编目数据集作为边界参考数据，检验冰川范围提取的准确性。气象数据选用离普若岗日冰原最近的班戈、安多、那曲、申扎和沱沱河 5 个气象站点的气温和降水量数据（https:/）。1.2 方法方法 1.2.1 影像预处理利用 GEE 中的筛选命令对 landsat 5/8 数据集进行时间、空间及云量的初步筛选，获取 19882020 年（除去 2012 年）59 月云量小于 50%27的研究区遥感影像集。基于 GEE 平台提取冰川面积的流程如图 2 所示。为解决仍有冰川区被云层遮盖的问题，采用像元级云掩膜算法，以无云像元替换有云像元，逐像素合成无云影像。为保证去云效果且高质量影像
13、能够被充分利用，每 3 a 合成一幅无云影像，即在研究时段内以相邻 3 a 为一期（其中 20132014 年为一期），每期数据选取 59 月进行影像合成和后续处理。图 3 为 20182020 年研究区影像云量去除前后对比，可见去云质量控制良好。1.2.2 冰川信息提取冰雪在绿波段具有较高的光谱反射值，而在短波红外波段光谱反射值很低，以这 2 个波段计算NDSI 可以有效区分冰雪与其他地物29。基于 GEE平台计算每期影像集的 NDSI，以各像元位置的最小 NDSI 值合成一幅新影像。最小 NDSI 合成法能够在非冰川区未被积雪覆盖时，排除该区域像素为冰川的可能性，从而在最大程度上减小季节性
14、积雪的影响27。NDSI 计算公式为NDSI=GreenSWIRGreen+SWIRGreenSWIR式中，为绿波段的反射率，在 Landsat 5 中对应波段为 B2，Landsat 8 中对应波段为 B3；为短波红外波段的反射率，在 Landsat 5 中对应波段 图 2 基于 GEE 平台的冰川提取技术流程Fig.2 Glacier extraction process based on GEE platform6 期田梦祺等：基于 Google Earth Engine 平台的青藏高原冰川变化研究以普若岗日冰原为例945 为 B5，Landsat 8 中对应波段为 B6。计算 NDSI
15、 前，在普若岗日冰川附近以 RGI6.0冰川编目数据为基准设置了一个 200 m 的缓冲区，初步排除研究区内积雪斑块等不确定因素的影响。根据多次实验和已有研究30，发现陆地卫星影像中冰雪的 NDSI 普遍高于 0.4，因此将 0.4 作为冰川提取的有效阈值。受积雪、阴影及冰碛等因素影响，影像中冰川内部会出现细小空隙或缺失像素，外围会存在离散斑块31。基于 GEE 平台，运用形态学开闭运算中的膨胀与腐蚀算法，可实现消除影像中细小空隙，平滑较大物体边界的效果，进而提高后续冰川边界提取及面积计算的准确性。1.2.3 冰川面积计算采用 GEE 中的碎斑去除算法，依据碎斑与目标冰川面积大小的差异（即周围
16、像素个数不一致），通过判断连通像素个数可实现不确定性斑块的去除，以减小像元误判和积雪斑块对提取精度的影响32。基于 ArcGIS 软件，并由栅格格式转换为矢量格式，进一步获得冰川面积。1.3 不确定性分析与精度评价不确定性分析与精度评价受季节性积雪、阴影以及影像分辨率的影响，利用遥感影像提取的冰川边界会存在误差，因此需要对提取结果进行精度评估，确保误差在合理范围内。研究采用 Bolch 等33提出的误差评估方法，对提取的冰川边界以遥感影像分辨率的 1/2(TM、OLI 影像为 15 m)作为缓冲区，将缓冲区边界内的冰川面积与原影像提取的冰川面积作差，以此差值与原始影像提取面积的比值作为评价指标。以19882020 年共 11 期冰川提取结果的不确定值为例，各期不确定值均在1.17%1.32%之间，处于合理的范围内19。2 结果分析 2.1 19882020 年普若岗日冰原面积变化年普若岗日冰原面积变化由 GEE 平台，可以快速获取 19882020 年普若岗日冰原的面积变化序列。与早期的研究结果对比，本文结果与前人研究结果28,34-36之间的绝对误差都在 5.6 km2以内，相对偏差