基于MODIS和CLCD数据的黄土高原土地利用变化检测及其对植被碳吸收模拟的影响

吴 巧 丽,张 鑫 阳,蒋 捷

(北京建筑大学测绘与城市空间信息学院,北京 100044;自然资源部城市空间信息重点实验室,北京建筑大学,北京 102616)

0 引言

为保护和修复生态系统,我国先后开展了三北防护林体系建设、天然林保护工程、退耕还林还草工程和沙漠化综合治理等重大生态工程[1]。黄土高原作为我国第一批退耕还林还草试点地区,近20年来为提高我国森林覆盖度作出了巨大贡献[2,3],1999—2019年黄土高原的植被覆盖率由31.6%提高到65%,有效遏制了该地区的水土流失问题[4],同时,实施退耕还林还草工程对于黄土高原地区生态系统碳汇能力提升具有重要作用[5,6]。目前,大量研究聚焦在土地利用/土地覆被变化(Land-Use and Land-Cover Change,LUCC)对黄土高原气候、生态环境和碳循环等的影响[7,8],缺乏对长时序土地利用土地覆被(Land Use and Land Cover,LULC)产品的应用准确性研究,制约着黄土高原退耕还林还草工程生态效益评估的准确性。

LULC产品是陆面生态过程模型进行植被—大气碳通量模拟的核心驱动数据,其产品精度直接影响区域尺度到全球尺度总初级生产力(Gross Primary Production,GPP)的模拟精度。目前,全球已经公开发布了多种不同时空分辨率的LULC产品,如MODIS LULC、CLCD LULC、GlobeLand和World Cover等[9,10]。然而,不同LULC产品的分类算法、样本选取和验证方法各异,导致产品应用于局部区域的精度存在差异,增加了区域尺度GPP估算的不确定性[11]。既有研究结果表明,MODIS LULC(MCD12Q1)数据在黄土高原地区的分类精度不足60%[9],导致基于MCD12Q1数据估算的GPP难以准确分析LUCC对黄土高原陆地生态系统碳循环的影响。

综上,本文引入时空分辨率和分类精度更高的CLCD LULC产品,对比分析CLCD LULC和MODIS LULC产品在表达2001—2019年黄土高原LUCC时空变化特征上的差异,并进一步基于长时序CLCD LULC产品、气象数据和多源遥感数据驱动GPP遥感机理过程模型,模拟黄土高原地区在退耕还林还草工程实施近20年来植被GPP的时空变化特征,定量分析LULC产品不确定性对黄土高原GPP估算结果的差异,以期提高黄土高原地区的碳汇估算精度,为精确量化国家重大工程生态效益以及政府调整和执行生态修复政策提供科学依据。

1 研究方法与数据来源

1.1 研究区概况

黄土高原涉及陕西、山西、甘肃、宁夏、河南、内蒙古和青海共7个省区,面积约6.4×105km2(图1),地势西北高东南低,地形复杂、沟壑纵横,主要地貌类型有黄土丘陵、黄土山地、黄土台地和河谷平原4类[12];植被类型从东南到西北依次为森林植被带、森林草原植被带、典型草原植被带、荒漠草原植被带和草原荒漠带[13]。该地区属于温带大陆性季风气候,年平均温度介于4～14 ℃,年均降水量在200～750 mm之间,跨越了半湿润、半干旱以及干旱区3个区域[14]。黄土高原恶化的生态环境严重制约当地社会经济可持续发展。为减缓水土流失,提高植被覆盖度,黄土高原实施了坡改梯、淤地坝建设和退耕还林还草等一系列重大生态建设工程,LULC格局发生显著变化[15-17],缓解了水土流失严重的问题,为提升陆地生态系统固碳能力奠定了坚实基础[18]。

图1 黄土高原地貌空间分布

1.2 LUCC变化像元检测和FGM模拟方案设计

为精确估算退耕还林还草工程实施近20年黄土高原地区陆地生态系统碳吸收的时空变化特征,本文首先基于不同源LULC数据,对比分析2001—2019年黄土高原LUCC时空变化特征,将研究期内发生土地利用变化的像元作为主要研究对象[19],若19年土地利用类别均一致,则识别为“稳定像元”,反之则为“变化像元”,提取变化像元制作土地利用转移矩阵;其次,基于MCD12Q1数据和CLCD数据,结合多源遥感数据和气象数据驱动Farquhar GPP Model(FGM)模型,生产2001—2019年黄土高原地区500 m空间分辨率、月尺度的GPP遥感产品,以CLCD GPP产品分析近20年黄土高原地区GPP的时空变化特征;最后,以CLCD LULC产品估算的GPP为基准,定量分析MCD12Q1分类误差对GPP空间分布、各类植被GPP年际变化特征、GPP年总量和年际变化趋势的影响。具体流程如图2所示。

图2 技术路线

1.3 FGM模型及其驱动数据

FGM模型[20]基于Farquhar模型进行叶片尺度的光合作用速率模拟,利用二叶模型进行叶片到冠层尺度的光合作用速率转换,并结合阴阳叶LAI比例估算冠层尺度GPP。FGM模型解决了Fick定律、Farquhar模型和Ball-Woodrow-Berry气孔导度模型联立求解胞间CO2浓度时,复杂度高、计算效率低、参数化难度高的问题,在保留模型机理性强优势的同时,提高了模型运算效率,更适合大尺度GPP模拟。目前,FGM模型已成功应用于中国[20]和欧洲[21]的GPP估算,因此,本文选用FGM模型定量分析LUCC对黄土高原GPP估算结果的影响,模型详细参数及运算过程参见文献[20-22]。

FGM模型的驱动数据包括遥感数据、气象数据和大气CO2浓度观测数据(表1)。其中,遥感数据包括MODIS LULC[23]和CLCD LULC产品[24]、全球陆表特征参量(GLASS)的叶面积指数(LAI)数据[25]和下行短波辐射(DSR)数据[26];气象数据来自国家地球系统科学数据共享服务平台—黄土高原科学数据中心,是中国学者基于机器学习算法、中国气象站点观测数据和CRUNCEP气象再分析数据集生产的一套空间分辨率更高的气象参数产品,主要包括气温、相对湿度以及饱和水汽压差(VPD)[27-29];由于缺乏时空连续的CO2浓度产品,本文采用美国夏威夷岛MLO观测站高精度和长时间跨度的大气CO2浓度实测数据,该数据可以捕捉大气CO2浓度的平均长期变化趋势[30]。

表1 FGM模型的驱动数据集

由于多源遥感数据的时空分辨率、投影坐标系、空间覆盖范围和存储格式均存在差异,为获取配套的模型驱动参数,需进行数据预处理,使不同输入数据的时空分辨率和投影一致,为GPP模拟奠定数据基础。本文将遥感数据转为tiff格式后,对数据进行拼接和裁减,重投影到WGS84地理坐标系,采用最邻近插值法将数据分辨率重采样到500 m;气象数据采用最邻近插值法将空间分辨率重采样到500 m,由于气象数据时间分辨率为月尺度,本研究将LAI和DSR产品重采样到月尺度。最终,本文利用表1数据驱动FGM模型,估算出2001—2019年黄土高原地区500 m空间分辨率的GPP月尺度数据,并将GPP产品和LULC数据均转投影为Albert_Conic_Equal_Area等面积投影坐标系,避免数据投影和面积变形对GPP年总量估算精度的影响。

1.4 土地利用数据及重分类方法

1)MODIS LULC数据包括MOD12Q1(1 km空间分辨率)和MCD12Q1(500 m空间分辨率)两种数据,本文选择较高空间分辨率的MCD12Q1产品。根据监督决策树分类方法,MCD12Q1产品包含5种不同的土地覆盖分类方案,本文按照国际地圈生物圈(IGBP)计划定义分为17类[31],包括常绿针叶林(ENF)、常绿阔叶林(EBF)、落叶针叶林(DNF)、落叶阔叶林(DBF)、混交林(MF)、郁闭灌丛(CSH)、稀疏灌丛(OSH)、多树草原(WSA)、稀树草原(SAV)、草地(GRA)、湿地(WET)、农田(CRO)、农田与自然植被混合地类(CRV)、冰雪(PSI)、城市建成区(URB)、荒地(BRN)和水体(WAT)。

2)CLCD LULC数据是由武汉大学杨杰等基于Google Earth Engine平台1982—2021年多传感器观测数据,利用随机森林分类器得到的LULC分类结果,空间分辨率为30 m[24]。在黄土高原地区CLCD LULC产品土地利用分类精度高于MODIS LULC产品[32]。针对MODIS和CLCD LULC产品分类体系的异同,参照全国土地利用遥感分类系统标准[33],本文结合黄土高原的地形地貌特点,将两套LULC分类体系按一级大类重分类为草地、林地、耕地、灌丛和非植被5种地类[20](表2)。两种LULC统一采用Albert_Conic_Equal_Area等面积投影坐标系,黄土高原空间范围矢量数据来源于中国科学院地理科学与资源研究所(https://www.resdc.cn/,2023年3月9日获得)。

表2 两种土地利用与土地覆盖(LULC)产品重分类前后的类别对应关系

2 结果分析

2.1 2001—2019年黄土高原LUCC时空分布特征

基于2001—2019年MODIS MCD12Q1和CLCD LULC数据得到的黄土高原各地类面积存在较大差异(图3)。MCD12Q1数据(图3a)显示,2001—2019年黄土高原耕地、林地和灌丛面积分别增加20 688 km2、8 993 km2和2 269 km2,但草地面积减少量(约26 635 km2)超过耕地和灌丛增加的总面积,与退耕还林还草实际不符,可能是MODIS分类方法混淆了耕地和草地所致;CLCD LULC数据(图3b)显示,2001—2019年黄土高原林地(+14 025 km2)和草地(+7 262 km2)面积增加,耕地(-18 843 km2)、灌丛(-1 727 km2)和非植被(-717 km2)面积减少,更符合退耕还林还草实际。

图3 基于MCD12Q1和CLCD LULC数据的2001—2019年黄土高原各地类面积变化

基于2001—2019年MCD12Q1和CLCD LULC数据得到研究区各地类变化像元空间分布(图4),统计各地类转换占比(图5)可知:①MCD12Q1数据(图5a)显示,研究期间黄土高原变化像元地类主要为草地(占变化像元总面积的42.47%)、灌丛(占比16.63%)和耕地(占比13.86%)。其中,草地主要转为耕地(25.97%)和灌丛(12.07%),灌丛主要转为林地(6.69%)和耕地(6.56%),少部分转为草地(3.09%),耕地主要转为草地(6.84%)和灌丛(4.77%)。此外,还有约8.07%的变化像元发生在非植被区域,主要转为草地(7.18%);林地像元变化最少(1.73%),主要转为灌丛(1.47%)。②CLCD LULC数据(图5b)显示,黄土高原变化像元地类主要为耕地和草地,分别占变化像元总面积的31.82%和31.00%,耕地主要转为草地(24.32%)和非植被(5.07%);草地主要转为耕地(18.93%)和林地(7.79%);还有9.36%的变化像元发生在非植被区域,主要转为草地(8.59%);林地像元变化较小,仅占变化像元总面积的1.26%,主要转为耕地(0.81%)。

图4 基于MCD12Q1和CLCD LULC数据检测的2001—2019年黄土高原变化像元的空间分布

图5 基于MCD12Q1和CLCD LULC数据检测的2001—2019年黄土高原土地利用转移矩阵

为更好对比两种土地利用产品在黄土高原的分类精度,本文在甘肃省庆阳市(107°38′46″E,35°43′51″N)选取典型耕地和草地混合地貌区。通过两种土地利用产品与Google Earth高分辨率影像对比,发现MCD12Q1确实对耕地和草地存在分类混淆现象,未能有效捕捉到Google Earth影像上的林地扩张现象;而CLCD LULC产品分类更精细且与Google Earth影像更接近(图6)。

图6 2001年和2019年黄土高原地区MCD12Q1产品、CLCD LULC产品和Google Earth影像对比

进一步根据分类精度更高的CLCD LULC数据计算2001—2019年黄土高原7个省区各地类变化面积(图7)。可以看出,7个省区耕地面积均有所减少,陕西、甘肃和山西位居前三,耕地面积分别减少6 146 km2、4 650 km2和3 228 km2;7个省区的林地面积均有所增加,山西、陕西和甘肃位居前三,分别增加6 303 km2、4 143 km2和1 945 km2;草地面积增量较大的为内蒙古、甘肃、陕西和宁夏,分别增加5 879 km2,2 289 km2、2 001 km2和1 593 km2,而山西草地面积减少最多,为4 933 km2。由此说明,退耕还林还草主要发生在黄土高原中部,尤其在黄河中游段西侧。

图7 基于CLCD LULC的2001—2019年黄土高原7个省区各地类面积变化

2.2 基于FGM模型的黄土高原总初级生产力估算

由上述研究可知,基于MODIS MCD12Q1和CLCD LULC数据得到的各地类空间分布、年际变化特征均存在显著差异,利用长时序CLCD LULC数据能更有效地捕捉退耕还林还草工程实施近20年来黄土高原地区LUCC的时空分布特征。因此,本文选取CLCD LULC产品估算2001—2019年黄土高原GPP,并以CLCD估算的GPP为标准,定量分析MCD12Q1数据分类误差对黄土高原GPP估算结果的影响。

由图8a可知,2001—2019年黄土高原GPP年均值为5 210.2 Tg C a-1,呈东南高、西北低的空间分布特征,最高值约为1 850 Tg C a-1,陕西、山西和河南GPP年均值明显高于其他省区;由图8b可知,林地GPP(年均值为104.2 Tg C a-1)对黄土高原GPP贡献最大,其次为草地(年均值为85.1 Tg C a-1)和耕地(年均值为82.6 Tg C a-1),灌丛贡献较小(年均值为2.4 Tg C a-1)。从年际变化趋势看,黄土高原地区陆地生态系统的碳吸收能力显著增强(年际趋势为+6.4 Tg C a-2,P<0.001),其中草地(年际趋势为+2.6 Tg C a-2,P<0.001)和林地(年际趋势为+2.5 Tg C a-2,P<0.001)贡献最大。

图8 基于CLCD数据和FGM模型估算的2001—2019年黄土高原GPP空间分布和各类型植被GPP年总量的年际变化趋势

2.3 LUCC不确定性对GPP估算结果的影响

图9a显示,基于MCD12Q1估算得到2001—2019年的GPP存在高值高估、低值低估的问题,即黄土高原南部地区(包括青海、甘肃、陕西省南部、山西、河南和宁夏南部)的MCD12Q1 GPP相对CLCD GPP偏大,北部地区(包括内蒙古、陕西北部和宁夏北部)MCD12Q1 GPP相对CLCD GPP偏小;估算结果显示,MCD12Q1 GPP相对CLCD GPP在林地存在明显低估现象(-38.3 Tg C a-1),在灌丛(+33.6 Tg C a-1)、草地(+28.9 Tg C a-1)和耕地(+2.3 Tg C a-1)存在高估现象;相对CLCD GPP估算结果,基于长时序MCD12Q1数据驱动FGM模型会给黄土高原地区GPP年总量模拟值带来9.6%的高估;从图9b看,2001—2019年MCD12Q1 GPP年总量约为5 712.8 Tg C a-1,年际增长趋势显著(+7.85 Tg C a-2,P<0.001),相比CLCD估算的GPP年际增长趋势产生22.7%的高估。

3 结论与讨论

本文基于MODIS MCD12Q1和CLCD两类典型土地利用数据,定量分析了2001—2019年黄土高原LUCC时空变化特征,揭示了MCD12Q1土地利用分类精度不足对陆地生态系统碳汇估算精度的影响,主要研究结论如下:

1)由于数据来源、分类标准和空间分辨率不同,MODIS MCD12Q1和CLCD两种LULC数据集显示的黄土高原LUCC的时空变化特征存在较大差异。黄土高原是退耕还林还草等一系列生态修复工程实施的重点区域,CLCD LULC数据空间分辨率更高,分类更精细,精度更高[24,32],该数据显示近20年来黄土高原耕地面积减少,林地和草地面积增加。在气候变化背景下,黄土高原各类植被GPP呈现显著增加趋势,证明退耕还林还草工程产生了重要的生态效益。值得注意的是,虽然MCD12Q1数据被广泛应用于全球GPP产品的生产,但其局地分类精度不一定高。例如,虽然MCD12Q1数据也显示出黄土高原林地面积呈逐年上升趋势,但表现出耕地面积增加和草地面积减少的特征,这与退耕还草的实际和谷歌高分影像目视解译结果均不相符,证明长时序MCD12Q1数据对黄土高原地区耕转林和耕转草解译精度较低,无法精确量化退耕还林还草工程的生态效益。

2)长时序CLCD LULC数据表明,2001—2019年黄土高原7个省区的土地利用变化规律基本一致,其中,山西和陕西退耕还林效果显著,而内蒙古以退耕还草为主,这与已有研究发现的草地退化规律一致[34,35],草地有可能进一步转化为林地[36]。综上,耕草转换是黄土高原的主要LUCC变化特征,部分地区存在耕转草、草转林的自然演替现象。

3)长时序CLCD数据能较准确捕捉黄土高原地区退耕还林还草现象,2001—2019年黄土高原GPP呈显著增长趋势,其中,林地和草地GPP增长趋势最显著,表明退耕还林还草工程对陆地生态系统碳吸收的促进作用。但植被光合作用速率受气候变化、大气CO2浓度增加、植被变绿、土地利用变化等众多因素共同作用的影响,需进一步分析植被碳吸收增长的具体原因。

4)研究假设CLCD LULC数据分类准确,定量分析了MCD12Q1数据分类误差对黄土高原GPP估算精度的影响。对比CLCD GPP,MCD12Q1 GPP在黄土高原南部和北部地区分别存在高估和低估现象,基于MCD12Q1数据估算的GPP年总量均值高估9.6%,GPP年际变化趋势则高估近22.7%;从各地类看,MCD12Q1 GPP低估了林地的GPP而高估了耕地、草地和灌丛的GPP。

虽然CLCD LULC数据在中国的分类精度更高[9,32],但目前无法完全避免土地利用分类误差,这也给GPP估算带来了一定的不确定性,未来有必要从土地利用分类误差向GPP模拟误差传递的角度开展进一步研究。