1　引 言

土地覆盖、土地利用是自然和人类力量共同塑造的结果，也是人类活动对地表自然生态系统影响的直接表现形式（Foley 等，2005；Turner 等，2007； Song 等，2018）。土地覆盖和土地利用变化，与陆地表层物质循环和生命过程密切相关，影响着生物圈与大气交互作用、生物多样性和生态系统服务功能、地表辐射能量平衡、生物地球化学循环以及资源环境的可持续利用（Yang 等，2013；Qiu 等，2020；Sy和Quesada，2020）。作为地球系统模式的关键输入参数之一，土地覆盖、土地利用数据对于全球环境变化及其模拟研究具有重要价值（Feddema 等，2005；Beusch 等，2020）。土地覆盖、土地利用数据也是农林生产、土地政策、城市建设、抵御洪涝和火灾防范的基础数据，可服务于生态资源环境和城市发展的评估、管理和决策，促进可持续发展（刘纪远 等，2014；刘芳 等，2016）。

21世纪以来，随着社会经济持续快速发展以及工业化、城市化进程的加快，中国耕地经受巨大压力，影响到中国粮食安全；同时一系列环境破坏问题，也对中国自然生态系统的可持续性构成威胁（Bai 等，2018）。为了紧守耕地红线和生态红线，需要高精度的土地覆盖、土地利用数据来监控土地变化。同时，解决城市快速扩张引起的高质量农地和其他生态系统损失和环境污染问题，以及制定城市规划、评估建设效果等都需要在较高的分辨率水平上对土地覆盖、土地利用进行识别（Gong等，2019；Gong等，2020）。此外，退耕还林、保护耕地等土地政策和生态保护工程的目标制定与实施效果评估，均需要多年份高精度土地覆盖、土地利用数据作为定量化的依据（Mao等，2019；Peng等，2019）。为了满足上述需求，至少需要逐年的空间分辨率在30 m甚至米级的土地制图更新数据；除了对空间分辨率有较高要求外，农业、生物多样性保护和能量平衡计算对土地覆盖、土地利用信息的时间频率要求要高于月甚至需要逐日水平，对地表覆盖类型的要求到作物或其他植被的种类水平。

国家统计局的统计数据，通常来自地方的调查数据，无法体现细致的空间差异性，且来自不同区域的数据受主观性影响，准确性难以评估。遥感卫星影像由于具备长时间监测能力、较高空间分辨率等特点，已成为土地覆盖、土地利用变化监测的重要数据源之一。为了满足各种应用需求，人们基于遥感数据开展了土地覆盖、土地利用制图的大量努力。然而任何一种产品都无法满足上述应用对空间分辨率、时间频率、土地覆盖/土地利用类型和分类精度的要求。现有粗分辨率制图产品包括全球5 km地表覆盖制图产品GLASS-GLC（Liu等，2020）、欧洲航天局气候变化倡议（ESA-CCI）300 m全球土地覆盖数据（http：//maps.elie.ucl.ac.be/CCI/viewer/index.php［2020-12-12］）、MODIS 500 m全球土地覆盖产品（Friedl等，2010；Sulla-Menashe等，2019），1 km国际地圈生物圈计划和信息系统覆盖图（IGBP-DISCover）（Loveland等，2000）、1 km马里兰大学土地覆盖图（Hansen 等，2000）、1 km全球土地覆盖2000（GLC2000）图（Bartholome和Belward，2005）等。高分辨率产品，如30 m和10 m精细分辨率的全球土地覆盖图（FROM-GLC，FROM-GLC10）（Gong等，2013，Gong 等，2019）、30 m全球土地覆盖数据产品（GlobeLand30）（Chen 等，2015）、中国多期土地利用遥感监测数据集（CNLUCC）（Liu 等，2005；Zhang 等，2014；Xu 等，2018）、美国国家土地覆被数据库（NLCD）（Vogelmann 等，2001；Homer等，2004，2007，2015，，2020；Xian 等，2009；Yang等，2018）和美国土地变化监测、评估和预测产品（LCMAP）（Brown 等，2020；Pengra 等，2020）。在全类别土地制图产品中，部分类别动态变化显著、制图精度低，出现了一些针对单类的土地制图产品，如水体（Wood 等，2011；Pekel 等，2016；Ji 等，2018；Pickens 等，2020）、不透水面（Pesaresi 等，2013；He 等，2018；Liu，Hu 等，2018；Xu等，2018；Gong 等，2020；Liu等，2020），森林（Hansen等，2013；Shimada等，2014；Crowther等，2015；Hamilton和Casey，2016），湿地（Hu等，2017；Murray等，2019）等，但不能满足反映多种土地类型变化的需求。

总体来看，已有的土地覆盖产品最高的时间频率通常到年，更高的时间频率，如季节尺度土地覆盖数据，很少（付薇 等，2014）。时序制图一致性、可比性差是一大难点（Sulla-Menashe等，2019；Zhu等，2020），降低了时序制图产品的可信度，阻碍了时序制图产品的应用（Gómez等， 2016）。已有大量研究开发了变化检测方法应对一致性问题，如（Kennedy等，2010；Verbesselt等，2012；Brooks 等，2014；Zhu和Woodcock，2014；Bullock等，2020；Zhu等，2020），然而缺少全类别的一致性变化探测技术。另一方面，常用的传统遥感数据源时间频率不够密集，因而难以有效反映时间频率更高的地表覆盖动态（Woodcock等，2020）。数据融合是一种有效应对方式（Zhang 2010；Wulder等，2018），通过融合多种数据，提高数据的时空分辨率，然而当前缺少将数据融合实际应用于大范围土地制图的应用。

样本问题是当前土地覆盖制图，尤其是大范围、高精度、高频率时序制图的阻碍之一（刘涵等， 2017；Liu等，2019；Woodcock等，2020）。土地覆盖、土地利用复杂多样，样本，尤其是时序样本，采集耗费大量人力、物力，成本高、难度大（Foody，2010；Zhao等，2014；Liu等，2017；Pengra等，2020）。然而通常建立的样本集往往只服务于当前制图任务和需要，样本的积累和共享不足（Li等，2017）。利用当前样本库进行迁移和复用方面的尝试偏少。目前，仅有FROM-GLC10一套全球土地覆盖产品对样本进行了跨传感器和时空尺度的迁移复用（Gong等，2019）。有限样本稳定分类等样本迁移理论的提出，为样本迁移复用、跨传感器、跨尺度快速制图提供了一定理论依据（Gong等，2019；Huang等，2020）。

遥感大数据是土地时序制图面临的另一大瓶颈（Baumann 等，2016）。遥感卫星数据量大、异构而复杂，存储、处理、分析需要大量工作和成本（Le等，2016；Baumann 等，2018）。传统的本地处理（例如台式计算机）和数据分发方法（例如基于场景的文件下载）无法适应遥感大数据（Baumann 等，2016）。为应对这些挑战，降低数据准备和处理的复杂性，出现了以数据立方体DC（Data Cube）为代表的新的遥感数据处理范式（Giuliani等，2017）。数据立方体旨在通过解决体积，速度和多样性挑战，以分析就绪的形式提供对大型时空数据的访问（Killough 2018）。当前，基于开放数据立方体等标准（Open Data Cube -https：//www.opendatacube.org/［2020-12-12］，Earth Observation Data Cube-http：//eodatacube.eu［2020-12-12］，Earth System Data Cube-http：//earthsystemdatacube.net［2020-12-12］，Earth on Amazon Web Services-EAWS-https：//aws.amazon.com/earth/［2020-12-12］和Google Earth Engine （Gorelick，Hancher et al. 2017）-GEE，https：//earthengine.google.com［2020-12-12］），已开发了一些DC。例如Digital Earth Australia （Lewis等，2017）， Swiss Data Cube （Giuliani 等，2017） 和Africa Regional Data Cube （Killough 2019）。这些应用往往依靠高性能计算机和云计算环境为数据基础架构提供支持（Lewis 等，2016；Rizvi 等，2018）。然而这些数据立方体往往只是将已有遥感卫星数据进行基本预处理，切片形式存储组织，并未对数据的缺失观测和异常值进行进一步的时空补齐，以便于时间序列分析。

Landsat是当前时间序列分析和制图中使用最广泛的卫星数据源（Woodcoc等，2008；Wulder等，2019）。这与Landsat数据开放获取、易于访问，校准和预处理相对成熟等特点有很大关系（Roy等，2010；Loveland和Dwyer， 2012；Wulder等，2012；Dwyer等，2018）。Landsat卫星数据30 m的空间分辨率也相对较高，然而与高空间分辨率相伴的是低时间频率。此外，Landsat卫星经历多次更新换代，传感器之间差异大。同时Landsat 7具有明显的条带现象，时序Landsat数据的时间一致性和数据质量难以保证。恶劣的大气条件（如厚云层），也常常使数据不完整，降低了数据的可用性。为应对这些缺陷和问题，减缓异常值和缺失值的影响，重建数据是一种有效的应对策略。据重建缺失信息时补充信息的不同来源，重建算法通常可以分为基于空间的方法、基于光谱的方法、基于时间的方法以及混合方法（Shen等，2015）。基于空间的方法使用数据的其余部分填充缺少的数据区域，代表性方法如内插方法（Zhang等，2007）。由于缺少足够的信息，基于空间的方法通常无法重建大面积丢失值。这一缺陷可以通过从光谱影像中的大量冗余光谱信息中重建特定频段中的丢失数据来克服（Wang等， 2006；Rakwatin等，2008；Shen等，2010）。但是，基于光谱的方法前提是损坏的光谱数据既具有不完整的又具有完整的光谱带。由于云的移动性和传感器的扫描偏差，在同一地理区域且在不同时期获取的数据可以提供补充信息。对于基于时间的方法，如果时间间隔太短，则两个连续数据集中的云将大部分重叠，并且时间相关性可能无效。但是，如果时间间隔太长，则土地覆盖范围可能会变化很大，从而破坏了相关性。常见有时间替代法（Melgani 2006；Chen等，2011；Zeng等，2013；Cheng等，2014），时间滤波法 （Roerink等，2000；Chen等，2004；Julien和Sobrino， 2010）。将这些方法相结合，利用在空间，光谱和时间域中的相关性可以取得更好效果，如（Benabdelkader和Melgani， 2008；Zeng等，2013；Cheng等，2014）。除重建方法外，已提出将MODIS与Landsat数据融合的多种算法（Gao等， 2006；Zhu等，2010），在填补缺失值的同时，还可以提高数据的时间频率。然而已有研究往往停留在算法研究领域，缺少大规模制图生产和使用的实验研究。

此外，已有的土地制图数据，常常针对土地覆盖或土地利用，无法满足其他应用需求。分类系统简单且灵活性差，无法根据不同应用需求，提供不同类别粒度的层次化产品（Yu等，2014；宫鹏 等，2016）。

考虑到已有数据和产品存在的缺陷和中国土地评估的重要性，为了把握和揭示中国土地覆盖、土地利用变化的时空特征，我们基于AWS云计算环境，生产了中国21世纪以来逐年多季节、多分类系统的30 m空间分辨率的土地覆盖、土地利用数据产品。为了提高Landsat数据一致性、时间频率和数据质量，提出基于虚拟星座的时空融合和重建技术，首次建成了21世纪中国逐日30 m的无缝数据立方体SDC（Seamless Data Cube）。我们开发了基于变化检测和时空一致性的长时序高分辨率动态制图框架。为了满足不同的应用需求，该产品中逐年数据的分类系统有3种类型、两级结构，土地覆盖分类系统有FROM-GLC分类系统、FAO Land Cover Classification System（LCCS），土地利用分类系统来自CNLUCC分类系统。针对FROM-GLC土地覆盖分类系统，生产了逐季节产品。

2　数据与方法

如前所述，气候变化研究、地球系统模式发展、及农林管理、碳氮循环和水循环模拟和生物多样性保护及防灾减灾等迫切需要遥感制图创新。这既是世界性的需求也是科学发展的必然趋势。信息技术发展到今天，使人类知识数字化，使人类智慧逐步在虚拟空间得以积累和应用，给既模仿学习人工遥感制图的灵活性、又利用计算机处理超高维多样化的海量数据的准确性和高效性创造了条件。人工智能技术中计算机的人脸识别精度超过人类代表着机器智慧达到了新高度。在遥感制图领域，我们把既能像专业制图人员一样综合利用色彩、纹理、相对位置、图像要素上下文、和其他辅助地理知识，又能充分发挥计算机在制图人员难以做到的快速、大规模处理超过三维空间的海量数据方面的优势的计算技术称为智慧遥感制图iMap（intelligent mapping of land cover and land use ）。作为实现智慧遥感制图的初步尝试，本文提出一套新的智慧遥感制图框架（图 1），将在这一框架下研制的第一版中国土地覆盖和土地利用逐年逐季节数据集作为智慧遥感制图1.0版本。

10.11834/jrs.20210580.F001智慧遥感制图概念示意图Intelligent remote sensing mapping

新智慧遥感制图框架旨在以高性能云计算平台为支撑，在综合利用国内外多源遥感数据、定量遥感产品和社交媒体等非传统遥感数据的框架下，构建虚拟星座，结合大数据驱动的机器学习等人工智能方法，探索和解决无缝智能数据重建、面向用户需求的层次化分类系统设计、普适样本库构建及迁移、时空一致性变化检测等基础性遥感难题，生产一系列高精度、具有国际先进水平的多尺度长时序全球地表覆盖制图产品，并提供基于云计算的在线制图平台，为遥感、地球系统科学以及社会经济发展提供关键性基础支撑，实现多尺度、多要素、精准化的智慧遥感制图服务。

在智慧遥感制图理念指导下，本文以21世纪中国为研究对象，提出了基于虚拟星座的无缝数据立方体融合重建的技术方案，和长时序高分辨率一致性动态制图框架（图 2）。依据ARD数据标准对Landsat和MODIS影像进行一致性预处理，组建虚拟星座，构成原始的数据立方体RDC（Raw Data Cube）。提出时空多源引导对抗修复网络等模型，对RDC进行影像质量修复与补全、时空融合和数据平滑，构建逐日30 m无缝数据立方体SDC（Seamless Data Cube）。基于SDC，进一步计算特征立方体（feature cube）。RDC、SDC和feature cube，共同构成多种层次的（multi-level）ARD产品，用户可根据需求直接应用。在制图方面，设计自适应用户需求的多种、多层次的土地覆盖、土地利用分类体系，以有限样本稳定分类理论为指导跨时间迁移使用普适样本库，基于自动机器学习进行模型调优和集成，然后以变化检测保证制图的时空一致性，最终构成长时序高分辨率动态制图框架。

10.11834/jrs.20210580.F002中国无缝数据立方体和地表覆盖制图流程The flowchart of seamless data cube reconstruction and land cover mapping for China

本文基于亚马逊云计算AWS（Amazon web Services）平台，搭建了在线实时、自动化、无服务器、端到端的遥感大数据生产链和并行制图系统，以完成全部生产处理流程。AWS公开数据集（AWS Open Data）提供了所需的Landsat、MODIS等遥感数据源。AWS高性能、高弹性、可扩展的分布式算力，满足了快速处理PB级数据的计算资源需求。借助AWS Step Functions，将AWS Lambda、Amazon ECS（Amazon Elastic Container Service ） 等云服务组合，搭建了自动化、无服务器、端到端的流式生产链。AWS提供了完善的人工智能和机器学习服务，采用Amazon SageMaker工具，高效设计和训练了影像修复模型和制图分类模型，并借助AutoML对模型结构和参数进行深度调优，最终完成分布式高性能推理。

2.1　基于虚拟星座的无缝数据立方体融合重建2.1.1　数据预处理

本文使用了来自亚马逊公开数据集（AWS open data）和谷歌云存储（Google cloud storage）的2000年—2019年所有可用的Landsat L1T （Level 1 terrain-corrected）影像。中国范围内一共包括579个图幅（Path/Row）。采用与Landsat分析就绪数据ARD（Analysis Ready Data）相同的范式对所有Landsat数据进行大气校正、生成表面反射率（Masek等，2006；Vermote等， 2016）、云及云阴影检测与掩膜（Zhu 和 Woodcock， 2014；Qiu 等，2019）等处理流程。在ARD处理范式外，我们还对影像进行了双向反射分布函数BRDF（Bidirectional Reflectance Distribution Function ）校正（Schaaf等，2002），生成BRDF校正反射率NBAR（Nadir BRDF-Adjusted Reflectance），以及对山区影像进行了地形校正（Soenenaf等，2005）。考虑到传感器之间的差异，我们对OLI和TM传感器进行一致性校正（Royaf等，2016）。所有数据组织和存储成空间规则、有时间标记、波段集中的网格（tile），得到原始的数据立方体RDC（Raw Data Cube）。Data cube组织范式优化了遥感数据的访问性能，具有时空高互操作性；时空可变网格，能以多种网格尺度分发；ARD处理标准使所有数据可直接用于分析。

为提高数据的时间频率和空间完整性，本文采用了来自AWS open data和美国地质调查局（USGS）同期的500 m 逐日MODIS MCD43A4 第六代NBAR数据产品与Landsat组成虚拟星座。用质量信息波段对MODIS进行预处理。所用波段如表1所示。按照同样范式，生成了MODIS raw data cube。

10.11834/jrs.20210580.T001

Landsat与MODIS对应波段

The corresponding bands of Landsat and MODIS
波段Landsat NBARMODIS NBAR
Landsat 5 TMLandsat 7 ETM+Landsat 8 OLIMCD43A4
Blue0.450—0.5200.450—0.5200.450—0.5100.459—0.479
Green0.520—0.6000.520—0.6000.530—0.5900.545—0.565
Red0.630—0.6900.630—0.6900.640—0.6700.620—0.670
NIR0.760—0.9000.770—0.9000.850—0.8800.841—0.876
SWIR 11.550—1.7501.550—1.7501.570—1.6501.628—1.652
SWIR 22.080—2.3502.090—2.3502.110—2.2902.105—2.155
μm2.1.2　影像质量修复与补全

MODIS数据观测频率高，互补信息相对多，采用经典的算法流程对MODIS影像进行修复。首先，本文利用邻近相似像元插补算法NSPI（Neighborhood Similar Pixel Interpolator）方法（Chenaf等， 2011）对MODIS缺失进行填补。NSPI方法原理简单，同时又较高的计算效率（Zhuaf等，2012）。NSPI方法基于同类地物光谱信息的空间连续性特征，以及这种特征在时间尺度上的不变性的合理假设，以邻近时间影像作为参考，使用加权线性模型预测目标像素（xw/2,yw/2,tP）

Fxw/2,yw/2,tp,B=?T1×∑i=1NWi1×Fxi,yi,tp,B+T2×Fxw/2,yw/2,tr,B+∑i=1NWi2×Fxi,yi,tp,B-Fxi,yi,tr,B

式中，F是反射率，B是波段，tp、tr分别是预测和参考时间，N是相似像素数量，w是窗口大小。?xw/2,yw/2?为预测中心像素，（xi,yi）是第i个相似像素位置。Wi1、Wi2是分别基于空间、光谱距离计算的相似像素的权重，T1、T2是将基于光谱—空间、光谱—时间的两种预测结合起来的权重。

进一步采用加权惠特克平滑算法（weighted Whittaker smoothing）（Eilers， 2003； Kong 等，2019） 平滑MODIS时间序列获得MODIS 重建的数据立方体（reconstructed data cube）。相比其他降噪和平滑算法，这种算法能够很好地平衡保真度和粗糙度。其本质是一种惩罚最小二乘算法，可最小化拟合误差并惩罚平滑曲线的粗糙度，目标是找到以下Q最小的最优平滑时间序列z：

Q=Swt+λRwt=FC-zTWwtFC-z+λDwtz2z=Wwt+λDTD/WwtFC

式中，Swt是保真度，Rwt是粗糙度，FC是原始时间序列，λ是粗糙度参数，Dwt是差分矩阵，Wwt是权重的对角矩阵，由对应点的质量决定。

Landsat数据纹理信息相对丰富，但观测频率低，影像修复相对MODIS要求高。为此，提出时空多源引导对抗修复网络（图 3），修复Landsat观测的条带、云和阴影等缺失和异常值。生成网络结构改进自（Lanarasaf 等，2018）提出的用于超分辨率Sentinel-2影像的DSen2，派生自EDSR（Lim等，2017），骨干是ResNet（He等， 2016）。主要包括卷积层、ReLU和过连接。采用长、加性的跳过连接（skip connection），将输入待修复影像与模型输出相加，这意味着网络学习预测的是对输入待修复影像的残差校正图。网络主体由D（网络深度）个ResBlock组成，每个ResBlock又由两个卷积层、一个ReLU层、一个残差缩放层（Res. Scal.）以及一个残差连接的加性层构成。网络输入包括待修复影像（当日Landsat影像）、缺失掩膜（当日Landsat影像条带、云和阴影掩膜）、上采样（upsample）的多源引导影像（当日MODIS影像）、季节时相剖面（待修复日期所在季节所有的影像分位数、物候等特征）。网络最终输出是当日完整的修复影像。判别网络部分，输入是修复影像与真实影像，主要包括卷积层、ReLU、全连接层和Sigmoid函数。为训练网络，构建了遵循全局随机分布的无云样本数据集，并在Landsat影像上模拟条带、云和阴影等噪声。假设y?是修复影像，y是真实影像，网络的损失函数?如下：

?=?MSE+α?PER+β?ADV

式中，?MSE=∑iyi?-yi1，是基于影像L1距离的像素级MSE损失，?PER=∑i?yi?-?yi22是测量特征空间中的距离的感知损失，??来自ImageNet数据集上预训练的VGG网络（Simonyan 和Zisserman 2014；Ren 等，2016），?ADV=∑ilog?(1-DI(yi?))则是对抗学习模块的对抗损失，有利于提高图像视觉效果（Ma等，2020）（He等，2017），DI表示判别网络。

10.11834/jrs.20210580.F003基于时空多源引导对抗修复网络的影像修复框架The framework for image reconstruction based on spatio-temporal multi-source guided adversarialreconstruction network

值得指出的是，所提出的时空多源引导对抗修复网络，可以结合多源的遥感数据进行互补修复引导，包括MODIS、SAR等。通过结合时空信息、对抗学习、残差学习等范式，能够高精度逼真修复影像任何条带、云或阴影导致的缺失。

2.1.3　时空融合

基于改进的时空自适应反射率融合算法ESTARFM（Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model）（Zhu等，2010），融合有观测日期重建的Landsat和MODIS数据。相比STARFM算法，ESTARFM通过考虑混合像元效应，极大提高了在异质区域的反演精度（Zhu等，2010），具有更好的适用性。对于给定区域，算法假设在同一时间获取的来自不同卫星传感器的遥感数据具有可比性，并且经过校正后彼此相关。同时考虑到具有相似反射率变化的相邻相同类别像素，利用移动窗口方法来充分利用来自相邻像素的信息。对每天，目标像素 （xw/2,yw/2,tp）Landsat反射率Fxw/2,yw/2,tp,B是两个就近参考观测日期tr1、tr2分别得到的预测的加权：

Fxw/2,yw/2,tp,B=Ttr1×Ftr1xw/2,yw/2,tp,B+Ttr2×Ftr2xw/2,yw/2,tp,B

式中，Ttr1、Ttr2为时间权重，依据MODIS对应参考时间和预测时间的的反射率变化计算。Fkxw/2,yw/2,tp,B是基于参考时间tk的预测值，k=r1,?r2，

Fkxw/2,yw/2,tp,B=Fxw/2,yw/2,tk,B+∑i=1NWi×Vi×Cxi,yi,tp,B-Cxi,yi,tk,B

式中，Vi是第i个相似像素的转化系数，根据粗分辨率像素内相似像素处的Landsat和MODIS的线性回归估计得到。Wi是第i个相似像素的权重，由空间距离和光谱距离计算得到。

对于融合结果进一步经过weighted Whittaker smoothing处理，得到长时序高一致性的高空间分辨率、高时间频率（30 m、逐日）的无缝数据立方体SDC（Seamless Data Cube）。

2.1.4　特征立方体

基于SDC，本文进一步构建了特征立方体（feature cube），计算了归一化植被指数NDVI（Normalized Difference Vegetation Index）（Townshend和Justice，1986），改进的归一化水体指数MNDWI（Modified Normalized Difference Water Index）（Xu 2006），归一化建筑物指数NDBI（Normalized Difference Built-up Index）（Zha等，2003），和归一化燃烧比NBR（Normalized Burn Ratio）（Miller和Thode，2007）4个常用的代表性指数。取各个波段每年的0，25，50，75，100分位数作为年际物候特征。本文按春季（3月、4月、5月）、夏季（6月、7月、8月）、秋季（9月、10月、11月）、冬季（12月、1月、2月）的划分，计算了各个波段每个季节的分位数作为季节物候特征。除光谱特征外，还采用多套窗口大小计算了空间特征，如灰度共生矩阵GLCM（Gray Level Co-occurrence Matrix）、拓展形态学属性剖面（EMAPs）（Song等，2014）等。其他辅助特征包括地理特征，经纬度，以及地形特征海拔、坡度、坡向。其中海拔数据来自美国航天局的30 m海拔数据SRTM（Shuttle Radar Topography Mission）。

2.2　长时序高分辨率地表动态制图框架2.2.1　分类系统

采用的分类系统涵盖土地覆盖、土地利用，共3种类型。土地覆盖分类系统来自 FROM-GLC分类系统（Gong等，2013）和LCCS（Bartholome 和Belward，2005），土地利用分类系统来自中国土地利用制图（CNLUCC）产品（Xu等，2018）所用的土地利用分类系统。这些分类体系均有两级结构、超过20种类别。FROM-GLC分类系统有10个一级类型，29个二级类型（L等，2017）。LCCS分类系统有8个一级类型，33个二级类型（L等，2017）。CNLUCC分类系统有6个一级类型，25个二级类型（Xu等，2018）。

2.2.2　样本集和样本迁移

土地覆盖样本来自首套全球30 m分辨率多季节样本集（Li 等，2017），具有FROM-GLC和LCCS两种分类系统标签。该样本集是多名有丰富经验的影像解译者从2015年左右获取的Landsat影像中收集，具有严格的解译和质量检查方案。其训练样本具有代表性和典型性，验证样本具有广泛性和系统性（Gong等，2013；Zhao等，2014）。全球首套10 m土地覆盖产品FROM-GLC10（Gong等，2019）基于该样本集生产。中国范围内，训练样本有5500个样本位置，验证样本有2100个样本位置，样本空间和类别分布分别如图 4（a）图 4（b）所示（审图号：GS（2020）7429号）。

中国制图样本空间分布（土地覆盖： 不同颜色代表不同的FROM-GLC一级类，各类比例如饼图;土地覆盖： 不同颜色代表不同的CNLUCC一级类，各类比例如饼图）Spatial distribution of mapping samples（Land cover： different colors represent different FROM-GLC level 1 classes, the proportions of each type are shown in the pie chart； Land use： different colors represent different CNLUCC level 1 classes, the proportions of each class are shown in the pie chart）10.11834/jrs.20210580.F004土地覆盖训练样本Land cover training sample10.11834/jrs.20210580.F005土地覆盖验证样本Land cover validation sample10.11834/jrs.20210580.F006土地利用训练样本Land cover training sample10.11834/jrs.20210580.F007土地利用验证样本Land cover validation sample

土地利用样本从CNLUCC产品（Xu 等，2018）中采集，21世纪以来包括2000年、2005年、2010年、2015年4期。先前研究已经证明了使用现有制图数据可作为分类训练的有效来源（Zhang和Roy 2017；Song等，2018；Brown等，2020；Liu等，2020）。为了保证样本的可靠性和可信性，进行了筛选。首先，样本必须来自类别的同质区域，避免类别边缘，以减少可能的类别混淆（Brown等，2020）。其次，当一个像素周围的8个像素相同时，中心的单元格才被认为是更可靠的样本（Wang等，2019；Brown等，2020）。为了使类别比例符合实际情况，本文从这些位置中，进行了分层随机采样。最终一共保留了46万随机点做训练样本，11万随机点做验证样本，训练和验证样本独立且符合同分布。样本数量级大，且具有多期标签，空间分布如图 4（c）—图 4（d）。

对于土地覆盖分类系统，将2015年分类器，用于所有年份分类。而土地利用分类系统，对每年，采用最邻近的1期分类器（2000年、2005年、2010年、2015年共4期）。尽管样本的迁移使用会引入一定错误，但有限样本稳定分类理论（Gong等，2019；Huang等，2020）指出了样本迁移使用的有效性。该理论显示，使用经全球制图反复试验证明表现突出的随机森林算法，在采用训练样本量仅占整体的40％，或训练样本的错误比例达到20％时，分类精度的损失仅为1%。

2.2.4　分类与后处理

Landsat数据空间分辨率是30 m，对空间上下文信息的依赖相对不高，此处没有采用基于卷积网络的语义分割类算法进行分类。考虑到重建的SDC像素层面的质量较高且信息丰富，我们采用基于像素的分类算法，从而使得结果能够保持更多丰富的纹理细节。基于自动机器学习AutoML（Automated machine learning）的思想，采用AutoGluon工具（Erickson等，2020）训练分类器，自动完成模型选择、超参数选择和调整等过程。采用的分类器包括随机森林（Random Forest）、超随机树（Extremely Randomized Trees）、k最近邻（k-Nearest Neighbors）、LightGBM提升树（LightGBM boosted trees）（Ke等， 2017）、CatBoost提升树（CatBoost boosted trees）（Dorogush等，2018）、神经网络（neural networks），并通过模型集成，进一步了提高分类精度。为减少分类结果中噪声干扰，我们对分类结果进行了后处理。采用最大后验概率马尔科夫随机场MAP-MRF（Maximum a Posteriori Markov Random Fields）（Liu和Cai，2012； Wang等，2015）对逐年分类标签和概率进行改进，减少空间噪声。同时采用基于季节和趋势分解的断点识别算法BFAST（Breaks For Additive Seasonal and Trend）（Verbessel等，2010）进行变化检测，每一段时间段内，将更新为该段最高置信度的类别，提高时序一致性和可比性。

3　结 果3.1　首套21世纪中国30 m逐日无缝数据立方体

基于提出的SDC重建技术，我们生产了首套21世纪以来中国逐日30 m分辨率SDC。结果展示如图 5所示（审图号：GS（2020）7429号），其中图 5（a）是2019年中值合成的假彩色影像，图 5（b）—图 5（e）是分别从2019年4个季节中随机抽取的一天的中国影像进行镶嵌。整体来说，SDC空间覆盖完整，具有空间异质性且较为合理，能很好地反映地表的真实状况。同时，SDC通过提供逐日影像，可以充分展现任一地点随时间变化的物候等动态信息，首次实现了高分辨率无缝遥感。此外，在SDC光谱波段外，为每幅影像配置了详细的质量控制（QA）波段，记录影像处理过程的所有信息，包括原始质量信息，如是否为云、阴影、缺失值等，以及重建处理过程和可信度等，具有高透明度和高可操作性。

中国无缝数据立方体假彩色合成示例Example of false color mosaic for China seamless data cube10.11834/jrs.20210580.F0082019年合成2019 composite10.11834/jrs.20210580.F0092019-04-152019-04-1510.11834/jrs.20210580.F0102019-07-152019-07-1510.11834/jrs.20210580.F0112019-10-152019-10-1510.11834/jrs.20210580.F0122019-01-152019-01-15

为了检验重建模型的精度和SDC的质量，本文从中国随机选择了500个tile进行检验。参考（Yan 和 Roy，2020）的方法，本文选择有效值≥20%的用于重建，从剩下有效值<20%的数据中选择有效值≥2%的用于检验。因此测试数据是独立的。此外，对测试数据进行了进一步滤波。考虑可能残留的部分云或阴影（Egorov等，2019；Qiu等，2019）影响，本文对测试数据的云和阴影mask进行20像素的膨胀，并舍弃这部分观测值，以尽可能减小云或阴影的发生率，保证测试数据质量。如图 6所示，SDC重建预测的反射率值与实际的观测值能很好的贴合，其中SWIR1、SWIR2波段相关系数R高，分别为0.91、0.90；Blue波段相对较低，为0.84，可能与蓝光波段受大气影响严重、质量相对差有关。以RMSE来说，平均水平为0.025左右，显示预测与实际观测的差距很小，具有较高的精度。就波段来说，近红外、SWIR1、SWIR2的RMSE相对其他波段高，与这些波段反射率水平偏高有关。就地区来看，东南地区RMSE水平最低，东北、华中地区较高。总体来看，本文提出影像重建技术方案具有高可行性，生产的高分辨率、高时间频率的SDC具有高可靠性、高精度的特点。

无缝数据立方体验证结果Seamless data cube validation results10.11834/jrs.20210580.F013验证数据集预测与观测反射率的散点图（黑色虚线是1∶1参考线）Scatter plot of predicted versus observed reflectance for the validation dataset（The black dashed line is the 1∶1 reference line）10.11834/jrs.20210580.F014不同地区验证数据集上反射率的RMSERMSE of predicted reflectance on the validation dataset in different regions3.2　首套21世纪中国30 m逐年逐季节土地覆盖、土地利用制图集

基于提出的长时序高分辨率地表覆盖制图框架，生产了首套21世纪以来中国30 m分辨率逐年土地利用和逐季节土地覆盖图集。2019年FROM-GLC、LCCS和CNLUCC等3种分类体系下制图结果依次显示在图7（a）、图7（b）、图7（c）（审图号：GS（2020）7429号）。从图7中可以看出，3种分类体系下反映的土地覆盖土地利用空间格局基本一致。西部、西北部地区的沙漠和草地地区，土地覆盖与土地利用分类体系下的制图结果有较大的差别。在西藏、内蒙古中部地区。土地覆盖结果中分为裸地，土地利用结果中分为草地。这种差异来自于土地覆盖与土地利用分类体系中对裸地、草地的定义、判别标准不同有关。2019年中国季节土地覆盖制图结果如图 8所示（审图号：GS（2020）7429号）。可以看到在FROM-GLC和LCCS两种分类体系下，季节动态变化都是显著的。春夏秋冬四个季节制图结果相比较，夏季的土地格局与年度制图结果更为相似。这与夏季植被生长态势旺盛，与其他地类更容易区分有关。同时，夏季受季节性冰雪影响小，而秋冬季节土地制图中，冰雪表征多。

2019年中国土地覆盖和土地利用图Land cover and land use map of China in 201910.11834/jrs.20210580.F015 FROM-GLC FROM-GLC 10.11834/jrs.20210580.F016 LCCS LCCS10.11834/jrs.20210580.F017 CNLUCC分类体系 CNLUCC classification system2019年中国季节土地覆盖图Seasonal land cover map of China in 201910.11834/jrs.20210580.F018春季Spring10.11834/jrs.20210580.F019夏季Summer10.11834/jrs.20210580.F020秋季Autumn10.11834/jrs.20210580.F021冬季Winter10.11834/jrs.20210580.F022春季Spring10.11834/jrs.20210580.F023夏季Summer10.11834/jrs.20210580.F024秋季Autumn10.11834/jrs.20210580.F025冬季Winter

局部区域结果显示了我们的制图结果可有效而细致反映年际土地变化动态。图 9显示的分别是重庆、西安的结果。重庆位于多云且复杂的山区，原始影像质量差，制图难度大。图 9显示我们的制图在重庆取得了准确而一致的结果，且对城市扩张、水体等有良好的刻画，而ESA CCI的结果缺少细节且不够准确。在西安，制图结果显示了细致的地表细节，尤其在村庄的识别上相比ESA CCI具有更高的准确性。

10.11834/jrs.20210580.F026中国重庆(西南)和西安(华中)制图结果比较（本文的：2000年—2019年，30 m；ESACCI：1992年—2015年，300 m）Comparison of mapping results for Chongqing (southwest China) and Xi’an (central China)（this paper: 2000—2019, 30 m; ESACCI:1992—2015, 300 m）

利用制图结果，统计中国21世纪以来土地变化情况如图10所示。从曲线来看，本文的制图结果年际波动小，时间一致性强。根据统计结果，耕地面积在2000年—2019年间呈显著线性下降趋势（p<0.001），下降速率为7684±814 km2/a （95%置信区间）。2019年中国耕地面积约为1.67×106 km2，未低于中国的耕地红线（约1.20×106 km2）。森林在19年间从2.29×106 km2增加到2.31×106 km2，增加速率为1670±176 km2/a（95%置信区间）。森林面积的显著增加（p<0.001）可能反映了中国退耕还林等生态工程的效果。不透水面是显著变化的另一大土地覆盖类别（p<0.001），2000年—2019年间约增加了134755 km2，相比2000年增加了约139%。

10.11834/jrs.20210580.F02721世纪以来中国各类地表覆盖面积变化Changes in the land cover area in China since the21st century

总体而言，本文的制图产品在全国范围内的各类曲线与其他不同来源制图产品和统计数据较为一致（图 10），得到的农地面积与5年间隔的CNLUCC产品接近，而300 m的ESACCI和2000年、2010年两期的Globeland30产品农地面积估计相对高。从趋势来看，除ESACCI显示农地增加，其他的时序制图产品农地趋势均为减少。此外，对农地在2015年的估计与30 m的GFSAD和FROM-GLC2015的结果也很接近。来自FAO和中国统计年鉴的统计数字，相对我们的产品以及其他制图产品，不连续，且对农地面积估计偏低，反映了统计数字与遥感制图产品的差异。这可能一定程度上导致了中国对耕地红线的过分担忧。对于森林，本文数据与ESACCI和CNLUCC估计接近，Hansen 2000年和Globeland30 2000年、2010年两期的估计相对低。中国统计年鉴2007年前数据与估计接近并略高于本文结果，而FAO估计整体偏低，两种统计数据来源均显示森林面积增加，与本文产品增加趋势一致。裸地和草地两大类别，由于各项产品间定义差别大，本文的估计与其他制图产品估计有差异。裸地和草地与具有同样分类体系的FROM-GLC2015年产品在2015年的面积更为接近。另一大主要类别，城市类别，可以看到中国统计年鉴的统计数字远远高于我们的估计和其他制图产品。城市类别面积与其他各种来源制图产品估计接近，略低于CNLUCC和Globeland30，而略高于GAIA和ESACCI。其中CNLUCC和Globeland30城市面积略高，主要是因为采用了人工解译方案，以及其对该类别的定义倾向于居民区。而ESACCI是所有制图产品中对城市类别估计最少的，这可能与其分辨率粗，识别出来的城市相对更少有关。尽管不同统计数字调查来源和方法不同，不同制图产品的分类体系、制图方法、数据来源各不相同，但总体来看，我们的产品各类的面积和趋势都处于合理的范围内。

图 11显示的是2019年分类不确定性图（审图号：GS（2020）7429号），可以看出，置信度高的地区往往是同质区域，如西北裸地、华东农地地区。西南地区置信度相对较低，可能与云雾条件有关，此外，黄土高原地区不确定性也相对较高，可能与位于裸地、草地、林地等交错边界有关。从类别来看，水体、裸地、冰雪、林地是置信度高的几类，而灌木、湿地两类相对较低，也反映了这两类制图难度大的特点。

10.11834/jrs.20210580.F0282019年中国土地覆盖分类置信度图(箱型图是每个类的置信度汇总)Classification confidence map in 2019(The box plot is a summary of the confidence for each class)

使用验证样本集对制图结果进行精度评价。包括2015年土地覆盖验证数据，4期土地利用验证样本。以及随机选择并解译的500个逐年土地覆盖、土地利用验证样本。评价指标包括整体精度OA（Overall Accuracy），卡帕系数（Kappa coefficient），各类用户精度UA（user’s accuracy）和生产者精度PA（producer’s accuracy）等。本文使用土地覆盖验证样本集对土地覆盖制图结果进行评估。

表 2结果显示2015年FROM-GLC一级类土地覆盖制图结果整体精度为80.8%，Kappa系数为0.75。其中林地、裸地两类PA、UA均超过了80%。农地类别PA、UA分别达到78.6%和84.6%。不透水面类别PA达到了75.0%。灌木、湿地两类由于样本比例较小，精度相对较低。灌木一类主要跟林地、草地两类混分，这3类定义根据植被高度和盖度进行划分，通过光谱信息划分难度通常较高。最新的LCMAP产品将草地、灌木两类合为一类。苔原一类，样本中数量极少，导致分类精度低。严格说来，分布在北纬60°以上地区，中国的范围实际并不属于苔原分布区。FROM-GLC二级类整体精度为73.3%，Kappa为0.67，除稀有类别外，各类别精度令人满意。

10.11834/jrs.20210580.T002

2015年FROM-GLC一级分类体系验证结果

Validation results of the 2015 FROM-GLC level 1 classification system
类别102030405060708090100UA/%
103082616800051084.6
203445091900010087.7
30232528123110356068.0
40647222800005025.9
50121038011017.7
600100111000084.6
7000100000100.0
8016121000361063.2
90412301202614594.2
1000020010051565.2
PA/%78.681.478.735.450.047.80.075.089.875.080.8

注：PA：生产者精度，UA：用户精度。

2015年LCCS一级类、二级类分类结果OA为79.0%、70.4%，kappa为0.72、0.64。耕地的PA、UA分别为77.3%、80.5%。林地、非植被的UA分别达到89.5%、94.3%。非植被的二级类人工表面及相关区域、水体两类PA、UA分别达到86.00%和84.6%、71.67%和95.5%。其中灌木/草本的UA相对较低，为65.2%，这再次说明灌木和草地分类难度大，分类效果不理想。

对CNLUCC土地利用分类系统数据评估结果显示，一级类2000年、2005年、2010年、2015年的OA分别为76.8%、77.0%、76.8%、76.0%。其中城乡、工矿、居民用地一类的UA偏低，这与CNLUCC人工解译方式和对该类定义倾向于居民区有一定关系。几类与湿地相关的类别的UA相对较低，与湿地类别分类难度大有关。湿地地区光谱时空变异性强，地表状态变化大，易与水体、草地、灌木等多类混分。

为了更好地掌握逐年制图精度，本文利用逐年样本对制图结果进行了检验。图 12中显示的是制图结果的OA。总体来看，分类总体精度在年际间有所波动。就土地覆盖结果来说，由于训练样本来自2015年，因而分类精度总体来看2015年相对最高。From-GLC、LCCS一级类的精度为80.0%±0.55%、78.4%±0.37%，level-2精度为72.0%±0.94%、70.0%±0.47%。CNLUCC土地利用分类结果2015年OA为77.2%，所有年份平均为77.2%。逐年精度结果总体令人满意。此外，为了检验基于变化检测的时间一致性处理的效果，我们同时检验了FROM-GLC一级类未经过该处理时的精度，结果如图 12中红色线显示。除2003年、2011年外，经过时间一致性调整的制图结果的精度均高于未进行处理的结果。这揭示了变化检测的有效性。

10.11834/jrs.20210580.F029逐年不同分类系统下制图结果精度Annual accuracy of mapping results under differentclassification systems

2015年FROM-GLC一级类分类体系下季节土地覆盖分类结果，春夏秋冬的OA分别为68.7%、70.5%、65.9%、66.9%。夏季分类精度最高，这可能与夏季为生长季，植物状态更好分辨有关。总体来看，季节制图结果精度低于年际制图结果，这与年际制图每年可用观测量多于一个季节的观测量有关。观测量少，云、大气状况、噪声等因素相对影响大，将严重干扰分类结果。

4　结 论

本文提出了智慧遥感制图的概念和一套智慧遥感制图的框架，该框架从用户需求出发、问题驱动，能够极大地改善当前数据产品难以满足农林管理、国情监测、生态环境保护、防灾减灾、城市建设等用户多样化、高精度地表监测需求的现状。

在该框架的指导下，本文基于AWS高性能、高弹性、可扩展的分布式计算资源，搭建了在线实时、自动化、无服务器、端到端的遥感大数据生产链和并行制图系统，并生产了首套21世纪中国全境逐日SDC及逐年逐季节土地覆盖和土地利用制图产品。逐日SDC，填补了高时空分辨率遥感观测的空白，作为ARD，对于高精度定量遥感反演和制图具有较大价值。逐年逐季节地表制图产品，平均精度超过80%，具有高时空一致性。这两套产品能够以前所未有的高空间分辨率和时间频率反映地表格局与动态变化，为科学研究、国民经济、社会发展提供了重要的基础地理数据，同时证明了智慧制图框架的可行性和有效性。未来，我们将进一步完善和发展智慧遥感制图理论和框架，以开放和灵活的理念，为促进中国遥感数据的易用和好用提供新思路。