生态学报  2023, Vol. 43 Issue (4): 1526-1536

文章信息

李作成, 张路, 欧阳志云, 胡金明
LI Zuocheng, ZHANG Lu, OUYANG Zhiyun, HU Jinming
基于Google Earth Engine的青藏高原土地沙化模拟与动态评估
Land desertification simulation and dynamic assessment of Qinghai-Tibet Plateau based on Google Earth engine
生态学报. 2023, 43(4): 1526-1536
Acta Ecologica Sinica. 2023, 43(4): 1526-1536
http://dx.doi.org/10.5846/stxb202203090566

文章历史

收稿日期: 2022-03-09
网络出版日期: 2022-10-13
基于Google Earth Engine的青藏高原土地沙化模拟与动态评估
李作成1,2 , 张路2 , 欧阳志云2 , 胡金明1     
1. 云南大学国际河流与生态安全研究院, 昆明 650000;
2. 中国科学院生态环境研究中心, 北京 100000
摘要: 青藏高原广布的高寒草原和高寒草甸,在气候变暖以及超载放牧等人类活动的影响下面临着沙化的风险,并且高寒干旱的气候条件使得高原草地生态系统极为脆弱,发生沙化后修复难度极大,对高原土地沙化时空动态评估极为必要。但目前针对大范围进行长时间快速沙化模拟监测的方法存在不足,通过对比实验构建了基于植被状况-地表反射率-土壤湿度(NDVI-Albedo-Wet)三维数据空间的沙化遥感模拟监测模型,利用Google Earth Engine遥感大数据对模型进行了代码实现,最终分析了青藏高原土地沙化时空格局的动态变化。结果揭示了青藏高原土地沙化总体缓解但问题依然存在的现象,具体体现在:(1)沙化总面积逐年减少,平均每年减少的面积达到55.13万hm2,但沙化模拟数据显示有大量的未沙化土地向轻度沙化转移。(2)沙化程度总体上趋于逆转状态,但部分区域土地沙化问题不断出现,主要分布于"一江两河"、唐古拉山南麓、柴达木盆地周边以及昌都地区,在三江源和川西甘南地区也有零星分布。
关键词: 青藏高原    土地沙化    模拟监测    时空格局    Google Earth Engine    
Land desertification simulation and dynamic assessment of Qinghai-Tibet Plateau based on Google Earth engine
LI Zuocheng1,2 , ZHANG Lu2 , OUYANG Zhiyun2 , HU Jinming1     
1. Institute of International Rivers and Ecological Safety, Yunnan University, Kunming 650000, China;
2. Eco-Environmental Research Center, Chinese Academy of Sciences, Beijing 100000, China
Abstract: The widespread alpine grasslands and alpine meadows on the Qinghai-Tibet Plateau face the risk of desertification under the influence of climate warming and human activities such as overloaded grazing. Moreover, the alpine and arid climate conditions make the plateau grassland ecosystem extremely fragile and difficult to repair after desertification. It is very necessary to evaluate the spatio-temporal dynamics of land desertification on the Plateau. However, the current method of large-scale rapid desertification simulation monitoring for a long time is insufficient. In this paper, a spatial desertification remote sensing monitoring model based on the vegetation condition-surface reflectance-soil moisture (NDVI-Albedo-Wet) three-dimensional data is constructed through comparative experiments, and the Google Earth Engine remote sensing big data platform is used to implement the model through code. Finally, the dynamic assessment of the long-term spatial pattern of land desertification on the Qinghai-Tibet Plateau is realized. The results reveal the overall mitigation trend of land desertification on the Qinghai-Tibet Plateau, but problems still exist. Specifically, it reflected that: (1) the total area of desertification is decreasing year by year, with an average annual reduction of 551300 hectares, but the desertification simulation data of each phase show that a large number of non-desertification land is transferred to mild desertification. (2) The degree of desertification tends to reverse on the whole, but the problem of land desertification in some regions continues to appear, mainly distributed in "one river and two rivers", the southern foot of Tanggula Mountain, the surrounding area of Tsaidam Basin and Qamdo area. It is also sporadically distributed in Sanjiangyuan, western Sichuan and southern Gansu.
Key Words: Qinghai-Tibet Plateau    land desertification    simulation monitoring    spatio-temporal pattern    Google Earth Engine    

全球有41%的陆地区域属于干旱区, 而全球农业用地中约有45%位于干旱区, 这些区域广泛的分布于中纬度地区, 在人为气候和土地利用等的影响下面临着沙化风险[1]。其中, 青藏高原因为其独特的地理位置, 具有独特的高寒干旱气候, 生态系统极度脆弱, 发生沙化后修复难度极大。高寒草地作为青藏高原的主要生态系统类型, 在过去几十年, 随着气候变化与人类活动加剧, 青藏高原超过50%的草地发生了不同程度的退化[24]。而青藏高原的植被一旦被破坏并产生裸露沙地, 不仅会给区域周边的生态环境和居民生产生活带来严重的危害, 而且地面粉尘上升后, 极易远程传输, 从而影响到整个西太平洋地区[56]。而目前对青藏高原土地沙化的研究大多局限于某一区域或者某一时间段进行, 缺乏对高原整体, 长时间序列的相关研究。

目前针对沙化监测已有较多的研究方法[711], 主要分为传统的地面调查结合遥感目视解译、较为成熟的基于各种指标阈值划分、以及计算机自动分类法(监督分类与非监督分类)。其中, 基于NDVI-Albedo特征空间划分土地沙化状态的方法, 近年来在各个地区都有所应用[1216]。但在面对较大范围或者较为复杂的生态环境, 应用NDVI-Albedo特征空间模型方法时, 常常存在“干、湿边”和数据的运算量以及匹配度问题。

针“干、湿边”问题, 有学者从多方面考虑尝试改善, 如李超[16]基于NDVI-Albedo特征空间, 通过作Albedo与NDVI的比率对模型进行优化, 希望能够改善原模型在水域监测方面的缺点。但在对“干、湿边”进行改善时, 较多的局限于NDVI-Albedo特征空间本身所依靠的Albedo与NDVI指数上, 没有考虑指数本身的生物物理意义是否完全合理。土地沙化的本质是土壤微生态系统总体退化的结果, 具体体现在植被盖度减少、土壤水分疏干、以及土壤颗粒粗化与土质疏松化。而Albedo不仅受植被覆盖影响, 同时受土壤湿度、地表粗糙度等的影响, 仅以NDVI-Albedo特征空间本身所依靠的Albedo与NDVI指数反映地表沙化过程或有不足。本研究通过对土地沙化相关遥感参量进行梳理, 进行对比实验来构建三维数据空间[17]替代NDVI-Albedo特征空间, 从而提高沙化模拟监测的精度。针对遥感数据的运算量以及匹配度问题, 在较大的区域尺度开展沙化遥感监测时, 常面临大量的遥感数据下载、校正、拼接等预处理过程, 以及多源遥感数据的时空匹配处理, 时间和人力成本较高。本研究利用Google Earth Engine遥感大数据平台在遥感数据处理方面的便捷性[18], 对实验建立的三维数据空间模型进行代码实现, 快速完成青藏高原区域, 多时相沙化信息的提取, 希望可以对土地沙化监测提供新的方法思路, 并为青藏高原土地沙化防治提供科学依据。

1 研究方法

基于NDVI-Albedo特征空间荒漠化差值指数方法, 考虑到该方法存在的“干、湿边”以及生物物理意义的明确性问题, 本研究拟通过几组具有较为明确生物物理意义的遥感参数指标, 进行对比试验来构建三维数据空间来替代NDVI-Albedo特征空间, 希望可以对沙化监测的准确度和精度有所改善。并在Google Earth Engine遥感大数据平台上通过代码实现, 以期能够实现青藏高原土地沙化长时间大范围的模拟观测。具体的试验流程及相关参数指标见图 1表 1

图 1 沙化监测模型构建流程图 Fig. 1 Construction flow chart of desertification monitoring model

表 1 沙化监测相关指标 Table 1 Indicators related to desertification monitoring
指标
Indicators
指标意义
Indicator meaning
指标计算公式
Index calculation formula
归一化植被指数
Normalized Difference Vegetation Index (NDVI)
利用植物叶片组织, 对蓝光和红光具有强烈的吸收, 以及对绿光尤其是近红外有强烈反射的特点。可以反映地表植被状况。 Landsat 5/7: NDVI = (B4-B3) / (B4+B3)
地表反照率Albedo[1921] 反映地表对太阳短波辐射反射特性的物理参量, 主要受太阳高度角、下垫面状况、土壤湿度以及气象条件等因素的影响。可以在一定程度上反映地表土壤的性状。 Landsat 5/7: Albedo = 0.356×B1+0.130×B3+0.373×B4+0.085×B5+0.072×B7-0.0018
遥感缨帽变换湿度指数Wet[22] 遥感缨帽变换所获得的湿度分量已被广泛地应用在生态环境监测中, 反映了水体和土壤、植被的湿度。可以在一定程度上反映地表土壤湿度。 Landsat 5/7: Wet=0.0315×B1+0.2021×B2+0.3102×B3+ 0.1594×B4-0.6806×B5-0.6109×B7
归一化水分指数
Normalized Difference Water Index (NDWI)[2324]
主要利用了植物在近红外波段具有最高的反射率, 而在中红外波段由于植物叶子水分的吸收作用导致反射率降低的特点, 可以反映地表在覆被状态下的水分状态。 Landsat 5/7: NDWI = (B5-B4) / (B5+B4)
表土粒径指数
Topsoil Grain Size Index (TGSI)[25]
表土粒径指数(TGSI), 与表层土壤的机械组成有关, 可以在一定程度上反映土壤质地。 Landsat 5/7: TGSI = (B3 - B1) / (B3 + B2 + B1)
地表温度T[22](Kelvin) 地表温度反映了地表与大气相互作用和大气与地面之间的能量通量, 是地表物理过程的重要参数之一, 本文中以校正后的热红外辐射通量表征。可以在一定程度上反映土地表层的热量状况。 Landsat 5/7: T = Scale×B6+Offset
B1、B2、...、B7分别代表Landsat Level 2, Collection 2, Tier 1地面反照率影像的波段1至波段7;Scale代表波段缩放因子,Offset代表波段偏移量
1.1 模型构建

(1) 基于NDVI-Albedo特征空间的荒漠化差值指数(DDI)原模型

基于NDVI-Albedo特征空间的荒漠化差值指数(DDI)方法, 是通过构建NDVI-Albedo特征空间, 即作NDVI与Albedo之间线性回归并通过回归线来表征沙化变化趋势, 然后在回归线上作垂线对NDVI-Albedo特征空间进行划分, 将不同程度的沙化土地有效的划分出来。考虑到NDVI-Albedo特征空间中沙化点位的分布, 本项研究中通过在研究区内划定不同沙化程度的典型区, 并在典型区内随机选取100个点, 提取其Albedo与NDVI指数后, 去除数据误差较大的点位共计76个数据点, 构建NDVI-Albedo特征空间如图 2, 对应的荒漠化差值指数(DDI)模型如下:

图 2 NDVI-Albedo特征空间图 Fig. 2 NDVI-Albedo feature space map

(2) 添加土壤湿度指数

考虑到在基于NDVI和Albedo指标的二维模型方法存在的“干、湿边”问题, 分别添加NDWI和Wet两种土壤湿度相关的遥感参量, 希望可以改善由地表湿度引发的模型监测精度问题。在整个研究区随机选取2500个数据点, 提取所需数据通过构建三维数据空间如图 3

图 3 基于NDVI-Albedo特征空间及土壤湿度(NDWI和Wet) 的三维数据空间图 Fig. 3 3D data space map based on NDVI-Albedo feature space and soil moisture (NDWI and Wet)

发现通过添加NDWI和Wet两种土壤湿度相关的遥感参量, 均可以较好的拟合实现进而对沙化状况进行监测, 具体模型如下:

(3) 以土壤湿度相关遥感监测参量与表土粒径指数(TGSI)代替Albedo

考虑到在基于NDVI和Albedo指标的二维模型方法, 以Albedo来表征土壤状况可能存在不足, 分别以两种土壤湿度相关遥感监测参量(NDWI和Wet)与表土粒径指数(TGSI)代替Albedo构建三维数据空间模型, 希望可以对土地沙化进行更加合理完善的表征。在整个研究区随机选取2500个数据点, 提取所需数据通过构建三维数据空间如图 4

图 4 基于土壤湿度(NDWI和Wet)、表土粒径指数(TGSI) 及NDVI的三维数据空间图 Fig. 4 3D data space map based on soil moisture (NDWI and Wet), topsoil particle size index (TGSI) and NDVI

发现基于NDVI、土壤湿度和表土粒径指数(TGSI)的三维模型拟合效果较差。其具体模型如下:

(4) 以土壤湿度相关遥感监测参量与地表温度(T)代替Albedo

考虑到在基于NDVI和Albedo指标的二维模型方法, 以Albedo来表征土壤状况可能存在不足, 分别以两种土壤湿度相关遥感监测参量(NDWI和Wet)与地表温度(T)代替Albedo构建三维数据空间模型, 希望可以对土地沙化进行更加合理完善的表征。在整个研究区随机选取约2500个数据点, 提取所需数据通过构建三维数据空间如图 5

图 5 基于土壤湿度(NDWI和Wet)、地表温度(T) 及NDVI的三维数据空间图 Fig. 5 3D data space map based on soil moisture (NDWI and Wet), surface temperature (T) and NDVI

发现基于NDVI、土壤湿度与地表温度(T)的三维模型拟合效果较差。其具体模型如下:

通过三维数据空间构建与线性拟合发现, 以土壤湿度相关遥感监测参量与表土粒径指数(TGSI), 或者以土壤湿度相关遥感监测参量与地表温度(T)代替Albedo构建三维数据空间模型对土地沙化进行监测的方法, 由于拟合优度太差不可取。而对基于NDVI和Albedo指标的二维模型方法通过添加土壤湿度指数构成的三维数据空间模型, 对提取的土地沙化数据有较好的拟合效果, 可进一步对比其与原模型沙化监测的准确性。

1.2 模型验证

利用Google Earth Engine遥感大数据平台对基于NDVI-Albedo特征空间和通过模拟实验得出的基于NDVI-Albedo特征空间及土壤湿度(NDWI和Wet)的三维数据空间模型进行代码实现。计算在各方法下高原区2020年期的沙化区分指数, 结合实际调查点位数据(有效点位共计885个, 大多于2020年秋季以及2021年夏、秋季节调查, 点位分布及沙化信息示例如图 6), 对土地沙化进行不同程度分级, 并进行精度验证。

图 6 实地调查样点分布及典型样地遥感图像 Fig. 6 Distribution of sample sites in field investigation and remote sensing images of typical sample sites
1.2.1 沙化模拟结果分级

以实际调查点位的位置信息, 提取GEE平台模拟得到的沙化监测差值指数。然后根据调查点位的沙化等级, 对所提取的沙化指数进行统计, 最终对模拟结果进行划分。

1.2.2 模型模拟精度验证

为了验证模型的监测精度, 结合实际调查点位数据的分级和各模型模拟结果的分级结果, 建立误差矩阵如表 2表 3表 4, 并计算Kappa系数[26]对各模型进行验证。Kappa系数计算方法如下:

表 2 DDI (基于NDVI-Albedo) 分类误差矩阵 Table 2 DDI (based on NDVI-Albedo) classification error matrix
分类
Classification
轻度沙化
Mild desertification
中度沙化
Moderate desertification
重度沙化
Severe desertification
极重度沙化
Extremely severe desertification
制图精度
Drawing accuracy
轻度沙化 62.69% 19.40% 11.94% 5.97% 62.69%
中度沙化 23.36% 52.80% 17.29% 6.54% 52.80%
重度沙化 12.39% 9.40% 55.13% 23.08% 55.13%
极重度沙化 4.69% 7.03% 18.75% 69.53% 69.53%
用户精度User accuracy 33.07% 71.97% 65.15% 55.28%

表 3 DDI (基于NDVI-Albedo-NDWI) 分类误差矩阵 Table 3 DDI (based on NDVI-Albedo-NDWI) classification error matrix
分类
Classification
轻度沙化
Mild desertification
中度沙化
Moderate desertification
重度沙化
Severe desertification
极重度沙化
Extremely severe desertification
制图精度
Drawing accuracy
轻度沙化 74.63% 22.39% 2.99% 0.00% 74.63%
中度沙化 12.62% 77.57% 9.35% 0.47% 77.57%
重度沙化 8.12% 11.54% 71.79% 8.55% 71.79%
极重度沙化 2.34% 10.94% 16.41% 70.31% 70.31%
用户精度User accuracy 50.51% 74.77% 79.62% 81.08%

表 4 DDI (基于NDVI-Albedo-Wet) 分类误差矩阵 Table 4 DDI (based on NDVI-Albedo-Wet) classification error matrix
分类
Classification
轻度沙化
Mild desertification
中度沙化
Moderate desertification
重度沙化
Severe desertification
极重度沙化
Extremely severe desertification
制图精度
Drawing accuracy
轻度沙化 88.06% 8.96% 2.99% 0.00% 88.06%
中度沙化 13.08% 71.50% 14.95% 0.47% 71.50%
重度沙化 1.28% 3.85% 87.18% 7.69% 87.18%
极重度沙化 0.78% 0.00% 4.69% 94.53% 94.53%
用户精度User accuracy 64.84% 91.07% 83.61% 86.43%

其中:OA代表每一等级类别正确划分的样本数量之和, 与总体样本数量的比值, 可以在一定程度上反映总体的分类精度。OB首先计算每一等级类别种, 实际划分到该类别的样本数量, 与模拟监测分级划分到该类别的样本数量的乘积, 然后根据每一等级类别得到的结果加和, 最后除以总样本数的平方, 可以用于校正数据的“偏向”问题。

通过比较原模型与分别添加NDWI和Wet两种土壤湿度相关的遥感参量构建三维数据空间模型的分类误差矩阵可知, 添加土壤湿度指数后, 分类的用户精度和制图精度均有不同程度的提高。结合原模型与添加土壤湿度指数的三维数据模型的Kappa系数(分别为0.43、0.63、0.77), 综合考虑拟合优度与Kappa系数, 最终确认基于Albedo-NDVI与Wet湿度数据的三维数据空间建立的沙化监测模型具有更好的适用性。

2 结果与分析 2.1 青藏高原相关数据收集及处理

青藏高原大部分处于中国西南部, 大致涵盖了中国西部、西北部两大沙区。包括西藏自治区和青海省的全部、四川省西部、新疆维吾尔自治区南部, 以及甘肃、云南的一部分。基于GEE平台, 利用Landsat(5/7) Collection 2 Level 2 Tier 1数据集和青藏高原边界的矢量数据, 筛选高原内6—8月生长季的遥感数据, 通过代码对数据进行去云处理, 然后对模型各指标(NDVI、Albedo以及Wet等)和最终沙化区分指数(DDI), 按表 1中参数计算公式和对比实验得到的沙化区分指数模型表达式进行测算, 最后以1988—2020年每三年数据, 根据NDVI最大值融合为11期沙化模拟数据。

2.2 青藏高原土地沙化现状空间分布格局

统计每个沙化强度等级相应的比例和总沙化面积, 发现2019年期(2018—2020)青藏高原沙化的总面积69.09万km2, 其中极重度沙化面积12万km2, 重度沙化面积16.89万km2, 中度沙化面积20.65万km2, 轻度沙化面积19.56万km2, 分别占沙化总面积的17.4%、24.4%、29.9%和28.3%。根据沙化现状分布图(图 7), 发现青藏高原土地沙化强度空间分异特征显著, 整体呈现“东南-西北”分异格局, 其中极重度和重度沙化主要分布在青藏高原西北部(青海省柴达木盆地、西藏阿里地区及新疆南部地区), 甘肃省南部和西藏自治区西南边界也有零星分布, 多属中国西部的沙漠聚集地带及其临近区域, 植被覆盖度低、气候干燥且多风;中度沙化主要分布在西藏自治区的西北部、柴达木盆地周边以及新疆维吾尔自治区南部, 在西藏“一江两河”区域也有分布, 多位于极重度和重度沙化区域边缘;土地沙化轻度的区域主要分布在西藏“一江两河”区域、青海省青海湖以及三江源区域, 处在“东南-西北分异线”上, 同时属于未沙化区域与沙化区域交错区域;无土地沙化区域主要分布在青藏高原的东南部、青海省东部及南部大部分地区。

图 7 青藏高原土地沙化现状(2019) 空间分布格局 Fig. 7 Spatial distribution pattern of land desertification status in the Qinghai-Tibet Plateau (2019)
2.3 青藏高原土地沙化时序分析

利用11期沙化模拟数据, 对模拟数据进行沙化等级划分、数据统计以及制图分析可得到如图 8的1990—2020年青藏高原沙化区域沙化时序变化图。从时序图中可以发现, 1990—2020年30间青藏高原极重度沙化和重度沙化面积总体上趋于稳定减少, 中度沙化面积有波动减少的趋势,轻度沙化呈现增加趋势。总体沙化面积由1989年期的8563万hm2递减到2019年期的6909万hm2, 平均每年减少的面积达到55.13万hm2, 沙化总体处于逆转趋势。其中, 极重度沙化向重度沙化、重度沙化向中度沙化、中度沙化向轻度沙化、轻度沙化与未沙化之间相互转变是主要变化量, 其次为重度沙化向极重度沙化转化、中度沙化向重度沙化以及中度沙化向轻度沙化转变。总体而言, 青藏高原沙化状况在程度上趋缓, 沙化处于逆转趋势, 但值得注意的是, 1990—2020年11期数据显示, 每一期沙化数据均有较为大量的未沙化向轻度沙化转移, 具体的变化量可以参考表 5

图 8 1990—2020年青藏高原沙化区域沙化时序变化图 Fig. 8 Time series change of desertification in the desertification region of the Qinghai-Tibet Plateau from 1990 to 2020

表 5 青藏高原土地沙化面积统计/(万km2) Table 5 Desertification area statistics of Qinghai-Tibet Plateau
年份
Year
极重度沙化
Extremely severe desertification
重度沙化
Severe desertification
中度沙化
Moderate desertification
轻度沙化
Mild desertification
总计沙化面积
Total desertification area
1989 20.71 24.46 22.80 17.66 85.63
1992 18.61 22.12 22.34 19.27 82.33
1995 16.78 21.20 22.57 18.86 79.42
1998 16.81 20.44 24.17 20.19 81.61
2001 14.05 18.55 23.08 19.32 75.01
2004 13.98 18.95 23.38 19.71 76.03
2007 14.15 17.69 23.75 20.58 76.17
2010 11.04 17.63 22.10 20.50 71.27
2013 16.37 19.08 21.13 18.51 75.09
2016 14.16 17.64 21.94 19.94 73.67
2019 12.00 16.89 20.65 19.56 69.09
2.4 青藏高原土地沙化动态变化空间分布格局

利用GEE平台以11期数据的年份为自变量、沙化区分指数(DDI)为因变量, 线性回归得到沙化区分指数(DDI)在30年间变化的斜率图(图 9), 可以反映1990—2020年青藏高原沙化区域土地沙化动态变化空间分布格局。

图 9 1990—2020年青藏高原沙化区域土地沙化动态变化空间分布格局 Fig. 9 Spatial distribution pattern of dynamic changes of land desertification in the desertification region of the Qinghai-Tibet Plateau from 1990 to 2020

通过1990—2020年青藏高原沙化区域土地沙化动态变化空间分布格局图可以看出, 青藏高原的沙化趋势整体上趋于逆转, 局部有沙化加剧现象。其中, 沙化逆转的区域面积约有86.4万km2, 面积占比约为69.3%, 沙化逆转明显的地区主要位于青海湖周边、三江源、那曲地区北部以及新疆南部区域;沙化加剧的面积约有38.3万km2, 面积占比约为30.7%, 沙化加剧明显的地区主要位于柴达木腹地、阿里地区北部、日喀则和那曲地区南部区域以及青藏高原东南缘地区。

对沙化加剧的区域, 通过叠加30年来未沙化与轻度沙化之间相互转变的数据与1990—2020年青藏高原沙化区域土地沙化动态变化空间分布格, 进一步分析发现:未沙化和轻度沙化相互转变区域有76.8%的面积是处于沙化逆转状态, 说明大部分未沙化向轻度沙化的转变是临时性的转变, 其转变原因可能是人类活动干扰或者短期气候变化扰动的结果;另外有23.2%的区域处于沙化加剧状态, 主要分布于“一江两河”、唐古拉山南麓、柴达木盆地周边以及昌都地区, 在三江源和川西甘南地区也有零星分布, 可能是接下来沙化发生的重点地区, 而其转移状态产生的具体原因可以进一步探讨。

3 讨论与结论 3.1 讨论

本文构建了基于Albedo-NDVI与Wet湿度数据的三维数据空间沙化监测模型, 并通过GEE平台进行代码实现。该模型相对于传统的沙化实地调查和遥感影像目视解译方法[8], 基于遥感指标进行模拟监测, 可以更加省时省力地实现大区域和长时间尺度的沙化监测。相较于指标划分、综合评价及地物提取中常用的机器学习算法, 本模型选用的三项光谱指数分别代表反照率(Albedo), 地表植被状况(NDVI)和土壤湿度(Wet), 具有更为明确的生物物理学意义[1011]。与Albedo-NDVI特征空间法相比[1216], 对土壤性状的考虑更为完善, 并根据模拟精度对土壤光谱参数进行了筛选, 由“干、湿边”问题造成的监测误差也得到了缓解。在实现方法上, 本文基于GEE平台建立三维数据模型, 充分利用了遥感大数据平台的便捷性[17], 相较于以往利用GEE进行地物分类方面的研究[2729], 也在一定程度上拓展了GEE平台的应用范围。总体而言, 本文通过模拟实验构建三维数据空间沙化监测模型, 进而利用GEE平台进行代码实现, 这一方法可以快速的对大范围土地沙化进行长时间序列的模拟监测。

根据三维沙化模型的长时序模拟结果, 本研究对青藏高原土地沙化近30年的时空动态特征进行了分析, 发现高原沙化总体上处于缓解的趋势, 但问题依然存在。在时间尺度上, 青藏高原土地沙化面积总体上呈现减少趋势, 这与以往国家荒漠化与沙化监测结果以及大家在青藏高原重点区域[3032](如三江源、青海湖等)监测的趋势相同, 但在沙化程度的年际转变中, 有中度及以上程度沙化和未沙化的土地向轻度沙化转变的趋势;在空间转变趋势上, 沙化程度整体上呈现逆转的趋势, 但还是存在较大面积的沙化程度加剧区域。一方面, 沙化加剧区域的存在说明了高原土地沙化问题依然存在, 另一方面, 通过对沙化加剧区域的分析探讨发现, 部分沙化加剧区域的转变是临时性的, 其转变原因可能是人类活动干扰或者短期气候变化扰动的结果, 这在与青藏高原气候以及畜牧政策等相关研究中也有体现[3334]

3.2 结论

本文通过土地沙化发生机制的文献回顾, 构建了基于总亮度-植被状况-土壤湿度(Albedo-NDVI-Wet)的三维土地沙化遥感监测模型。并利用Google Earth Engine遥感大数据平台进行了代码实现。对青藏高原土地沙化进行了长时间序列(1990—2020年)遥感模拟监测, 结果揭示了青藏高原土地沙化在总体缓解的趋势下, 问题依然存在。具体体现在沙化面积总体减少和沙化程度总体趋缓的同时, 部分区域的土地沙化问题不断出现, 主要分布于“一江两河”、唐古拉山南麓、柴达木盆地周边以及昌都地区, 在三江源和川西甘南地区也有零星分布, 建议未来生态修复工程重点关注以上区域。

参考文献
[1]
Burrell A L, Evans J P, de Kauwe M G. Anthropogenic climate change has driven over 5 million km2 of drylands towards desertification. Nature Communications 11, 3853(2020).
[2]
Chen H, Zhu Q A, Peng C H, Wu N, Wang Y F, Fang X Q, Gao Y H, Zhu D, Yang G, Tian J Q, Kang X M, Piao S L, Ouyang H, Xiang W H, Luo Z B, Jiang H, Song X Z, Zhang Y, Yu G R, Zhao X Q, Gong P, Yao T D, Wu J H. The impacts of climate change and human activities on biogeochemical cycles on the Qinghai-Tibetan Plateau. Global Change Biology, 2013, 19(10): 2940-2955. DOI:10.1111/gcb.12277
[3]
张镱锂, 祁威, 周才平, 丁明军, 刘林山, 高俊刚, 摆万奇, 王兆锋, 郑度. 青藏高原高寒草地净初级生产力(NPP)时空分异. 地理学报, 2013, 68(9): 1197-1211.
[4]
陈槐, 鞠佩君, 张江, 王元云, 朱求安, 颜亮, 康晓明, 何奕忻, 曾源, 郝彦宾, 王艳芬. 青藏高原高寒草地生态系统变化的归因分析. 科学通报, 2020, 65(22): 2406-2418.
[5]
孙东怀, 鹿化煜. 晚新生代黄土高原风尘序列的粒度和沉积速率与中国北方大气环流演变. 第四纪研究, 2007, 27(2): 251-262. DOI:10.3321/j.issn:1001-7410.2007.02.010
[6]
方小敏, 韩永翔, 马金辉, 宋连春, 杨胜利, 张小曳. 青藏高原沙尘特征与高原黄土堆积: 以2003-03-04拉萨沙尘天气过程为例. 科学通报, 2004, 49(11): 1084-1090. DOI:10.3321/j.issn:0023-074X.2004.11.013
[7]
王晓慧. 沙化土地遥感监测机理和方法研究[D]. 北京: 中国林业科学研究院, 2007.
[8]
王荣宝, 张楠楠, 王倩倩. 沙化土地遥感监测方法研究及应用. 测绘与空间地理信息, 2016, 39(10): 163-166. DOI:10.3969/j.issn.1672-5867.2016.10.047
[9]
徐玲玲, 延昊, 钱拴. 基于MODIS-NDVI的2000-2018年中国北方土地沙化敏感性时空变化. 自然资源学报, 2020, 35(4): 925-936.
[10]
胡梦甜, 张慧, 高吉喜, 仇宽彪, 王延松, 吕久俊, 鞠昌华. 基于RWEQ模型修正的土地沙化敏感性评价. 水土保持研究, 2021, 28(1): 368-372.
[11]
霍艾迪, 张广军, 武苏里, 刘志丽. 国内外荒漠化动态监测与评价研究进展与存在问题. 干旱地区农业研究, 2007, 25(2): 206-211.
[12]
曾永年, 向南平, 冯兆东, 徐豁. Albedo-NDVI特征空间及沙漠化遥感监测指数研究. 地理科学, 2006, 26(1): 75-81.
[13]
任艳群, 刘海隆, 唐立新, 姜亮亮, 安小艳. 基于NDVI-Albedo特征空间的沙漠化动态变化研究——以准格尔盆地南缘为例. 水土保持通报, 2014, 34(2): 267-271, 325.
[14]
岳辉, 刘英. 基于NDVI-Albedo特征空间的陕西省干旱与荒漠化遥感监测. 西北林学院学报, 2019, 34(1): 198-205.
[15]
李宇君, 张磊. 基于沙地指数模型的沙地监测方法. 地球信息科学学报, 2021, 23(4): 680-691.
[16]
李超. 基于差值指数的毛乌素沙地荒漠化监测分析[D]. 抚州: 东华理工大学, 2019.
[17]
薛丽红. 三维空间点中基于最小二乘法的分段直线拟合方法. 齐齐哈尔大学学报: 自然科学版, 2015, 31(4): 84-85, 89.
[18]
郝斌飞, 韩旭军, 马明国, 刘一韬, 李世卫. Google Earth Engine在地球科学与环境科学中的应用研究进展. 遥感技术与应用, 2018, 33(4): 600-611.
[19]
肖登攀, 陶福禄, Moiwo Juana P. 全球变化下地表反照率研究进展. 地球科学进展, 2011, 26(11): 1217-1224.
[20]
Liang S L. Narrowband to broadband conversions of land surface albedo I: Algorithms. Remote Sensing of Environment, 2001, 76(2): 213-238.
[21]
Liang S L, Shuey C J, Russ A L, Fang H L, Chen M Z, Walthall C L, Daughtry C S T, Hunt R Jr. Narrowband to broadband conversions of land surface albedo: Ⅱ. Validation. Remote Sensing of Environment, 2003, 84(1): 25-41.
[22]
徐涵秋. 区域生态环境变化的遥感评价指数. 中国环境科学, 2013, 33(5): 889-897.
[23]
Gao B C. NDWI-A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sensing of Environment, 1996, 58(3): 257-266.
[24]
邢丽玮, 牛振国, 张海英. 不同湿度指数在湿地分类中的对比研究. 地理与地理信息科学, 2015, 31(6): 35-40, 46.
[25]
Liu Q S, Liu G H, Huang C. Monitoring desertification processes in Mongolian Plateau using MODIS tasseled cap transformation and TGSI time series. Journal of Arid Land, 2018, 10(1): 12-26.
[26]
许文宁, 王鹏新, 韩萍, 严泰来, 张树誉. Kappa系数在干旱预测模型精度评价中的应用——以关中平原的干旱预测为例. 自然灾害学报, 2011, 20(6): 81-86.
[27]
Pekel J F, Cottam A, Gorelick N, Belward A S. High-resolution mapping of global surface water and its long-term changes. Nature, 2016, 540(7633): 418-422.
[28]
郭永强, 王乃江, 褚晓升, 李成, 罗晓琦, 冯浩. 基于Google Earth Engine分析黄土高原植被覆盖变化及原因. 中国环境科学, 2019, 39(11): 4804-4811.
[29]
Shang Y Q, Zheng X Q, Han R Q, Liu W C, Xiao F. Long-term evaluation on urban intensive land use in five fast-growing cities of Northern China with GEE support. Scientific Reports 11, 2073, 4(2021).
[30]
屠志方, 李梦先, 孙涛. 第五次全国荒漠化和沙化监测结果及分析. 林业资源管理, 2016(1): 1-5.
[31]
胡光印, 董治宝, 逯军峰, 颜长珍. 黄河源区1975-2005年沙漠化时空演变及其成因分析. 中国沙漠, 2011, 31(5): 1079-1086.
[32]
王静洁, 蔡延玲. 青海湖及其周围地区沙化土地变化动态研究. 湖南林业科技, 2016, 43(6): 81-89.
[33]
赵海霞, 那日苏, 刘雅学, 张文静, 秦艳. 草原植被退化与气候波动跨学科研究评述. 科技创新导报, 2011, 8(4): 144-145.
[34]
冯云飞, 李猛, 李少伟, 邸迎伟, 沈振西, 张宪洲, 余成群, 严俊, 席永士, 武建双. 2010-2017年藏北高寒退化草地禁牧恢复效果评价. 草业科学, 2019, 36(4): 1148-1162, 921.