About Earth Timelapse

地球延时摄影数据集提供了 1984 年至 2022 年这 40 年间地球变化的直观记录。该数据集由 Google 根据 NASA/USGS Landsat 任务（Landsat 4、5、7、8 和 9）和 ESA 的 Copernicus Sentinel-2 星群获取的 PB 级机载观测数据合成，以年度镶嵌图的形式提供，用于直观解读。

年度镶嵌图是全球完整、填补了空白且经过水体遮盖处理的可视化底图，可为公开的 Google 地球延时摄影互动查看器提供支持。使用像素级时间线性回归插值在相邻的有效观测年份之间平滑地重建 1999 年之前的观测空白。开阔的海洋水域替换为样式化的 NOAA ETOPO1 阴影浮雕测深数据和全球水域遮盖数据（使用 Hansen GFC 和 MOD44W 数据集）。这些数据适用于：直观解读、教育故事讲述、基础地图绘制或自定义视频导出。

空间分辨率时代

地球延时摄影集涵盖了两个不同的空间分辨率时代：

  • 1984 年至 2014 年（30 米 / 像素）：根据 USGS/NASA Landsat 4、5、7 和 8 的观测结果合成。图片以 30 米分辨率在 Web 墨卡托投影 (EPSG:3857) 中进行网格化。
  • 2015 年至今（19.11 米 / 像素）：通过将 ESA Copernicus Sentinel-2A/2B 多光谱成像仪 (MSI) 观测数据（10 米 / 20 米原生分辨率）与 Landsat 8 和 9 的数据融合而成。图片以 Web 墨卡托投影 (EPSG:3857) 方式网格化，分辨率约为 19.11 米。

年度镶嵌图覆盖了全球所有陆地和沿海地区，纬度范围约为 -82.6°S 至 83.69°N。要整合 40 多年的地球观测数据，需要协调不同卫星星座在运行寿命、传感器模式、轨道几何形状和大气扰动方面的差异。

地球延时摄影流水线会从美国地质调查局/美国航空航天局 (USGS/NASA) Landsat 计划（Landsat 4、5、7、8 和 9）和欧洲航天局 (ESA) Copernicus Sentinel-2 卫星群中提取数百万个单独的场景。为了将原始观测数据和大气层顶部 (TOA) 观测数据转换为一致的年度基图，系统会使用以下方法处理每个场景：

  • 物候和季节性场景过滤，以优化植被生长高峰期并最大限度地减少短暂的季节性影响。
  • 缺少通道恢复和传感器异常缓解（例如，Landsat 8 TIRS 饱和、Landsat 7 SLC-off 制品）。
  • 多传感器云、雾霾和阴影质量评分，结合了 Landsat simpleCloudScore、启发式 HSV 空间云恢复和 Sentinel-2 Cloud Score+。
  • 全色归一化和结构性伪全色锐化，以提高表观分辨率。
  • 使用多尺度低通空间滤波消除场景边界，同时保留局部地表纹理的 MODIS BRDF 校准辐射归一化。
  • 注重节水的边缘锐化，避免海岸线漂白。
  • 时间性 16 位中值合成，用于滤除瞬时大气噪声、阴影和短暂的伪影。
  • 时间线性回归插值（适用于全球最终版资源），用于重建 1999 年之前的镶嵌影像中缺失的观测数据。
  • 全球海洋水深测量建模和水体掩膜，结合了 NOAA ETOPO1、Hansen 全球森林变化数据掩膜和 MODIS 水体掩膜。
  • 全局色彩平衡、极地冰白点保留、多尺度局部对比度增强 (LCE)，以提升视觉效果。

源数据和传感器规范

Earth Timelapse 系列利用了以下机载地球观测任务：

卫星和传感器规范

卫星任务 传感器 运营日期范围 空间分辨率（MS / Pan） 提取的光谱频段 Earth Engine 集合 ID 马赛克分辨率时代
Landsat 4 Thematic Mapper (TM) 1982-07-01 至 1993-12-14 30 米/不适用 蓝、绿、红、近红外 (NIR)、短波红外 1 (SWIR1)、热红外、短波红外 2 (SWIR2) LANDSAT/LT04/C02/T1 30 米（1984 年至 1993 年）
Landsat 5 Thematic Mapper (TM) 1984-03-01 至 2012-12-31 30 米/不适用 蓝、绿、红、近红外 (NIR)、短波红外 1 (SWIR1)、热红外、短波红外 2 (SWIR2) LANDSAT/LT05/C02/T1 30 米（1984 年至 2012 年）
Landsat 7 Enhanced Thematic Mapper Plus (ETM+) 1999-04-15 至 2013-08-31 30 米 / 15 米 蓝、绿、红、近红外 (NIR)、短波红外 1 (SWIR1)、热红外、短波红外 2 (SWIR2)、全色 (Pan) LANDSAT/LE07/C02/T1
LANDSAT/LE07/C02/T2		 30 米（1999 年至 2013 年）
Landsat 8 Operational Land Imager (OLI) / TIRS 2013-02-11 - 至今 30 米 / 15 米 近岸蓝色、蓝色、绿色、红色、近红外 (NIR)、短波红外 1 (SWIR1)、短波红外 2 (SWIR2)、全色、卷云、热红外 LANDSAT/LC08/C02/T1_RT_TOA
LANDSAT/LC08/C02/T2_TOA		 30 米（2013 年至 2014 年）
19.11 米（2015 年至今）
Landsat 9 Operational Land Imager 2 (OLI-2) / TIRS-2 2021-11-01 - 至今 30 米 / 15 米 近岸、蓝色、绿色、红色、近红外 (NIR)、短波红外 1 (SWIR1)、短波红外 2 (SWIR2)、全色、卷云、热红外 LANDSAT/LC09/C02/T1_TOA 19.11 m（2021 年及更高版本）
Sentinel-2A / 2B 多光谱成像仪 (MSI) 2015-06-23 - 至今 10 米 / 20 米 / 60 米 蓝色、绿色、红色、RedEdge（1-4）、近红外、短波红外 1、短波红外 2、水汽、卷云 COPERNICUS/S2_HARMONIZED 19.11 米（2015 年及之后）

表 1. 核心光学传感器系统已集成到 Earth Timelapse 年度镶嵌图制作流水线中。

辅助数据集和参考数据集

  1. MODIS 全球 BRDF / 地表反射率基准：
    • 预先计算的多年地表反射率合成数据源自 Terra 和 Aqua MODIS（MODIS/006/MCD43A4 和 MOD09GA/MYD09GA），分辨率为 500 米，可作为行星辐射参考目标，用于大规模光照和雾霾协调。
  2. NOAA ETOPO1 全球地形：
    • 利用 1 角分的全球地形和水深数据，在全球水体掩膜产品中生成逼真的阴影浮雕海底地形和深度轮廓。
  3. Hansen 全球森林变化 (UMD/hansen/global_forest_change_2015)：
    • datamask 图层（区分陆地、永久性水体和沿海地区）是主要的陆地和水体分离边界。
  4. MODIS 全球陆地/水体掩码 (MODIS/MOD44W/MOD44W_005_2000_02_24)：
    • 与 Hansen GFC 结合使用，以隔离陆地表面并防止海岸水掩膜溢出。
  5. Sentinel-2 Cloud Score+ (GOOGLE/CLOUD_SCORE_PLUS/V1/S2_HARMONIZED)：
    • 像素级质量评估，可为 Sentinel-2 MSI 数据提供 cs（晴空置信度）和 cs_cdf（累积分布晴空概率）指标。
  6. 全球云气候统计信息：
    • 根据经验得出的第 25 百分位和最低云得分基准，可用于调整持续多云的热带地区的云掩盖阈值。

空间和时间过滤

全球年度镶嵌图可平衡以下需求：在多云的热带地区获取足够多的有效观测结果，同时防止高纬度地区出现冬季积雪污染和极端的太阳天顶角。

纬度和物候数据选取

  • 北半球极地 / 高纬度地区（>60°N 至 83.69°N）：
    • 按一年中的第几天 (DOY) 直接过滤：观测结果仅限于 DOY 150 至 270（5 月下旬至 9 月）。这以植被生长旺季（最大绿度）为目标，并最大限度地减少了因太阳高度角较低（太阳高度角 ≤ 0°）而造成的季节性积雪、冰盖和长地形阴影。
  • 温带和南半球地区（南纬 57° 至北纬 60°）：
    • 利用整个日历年（1 月 1 日至 12 月 31 日）最大限度地提高场景可用性。

处理方法：已遮盖的年度镶嵌

经过掩码处理的年度镶嵌流水线会生成未插值的全球合成影像。

缺少渠道恢复和传感器准备

在进行大气校正和质量评分之前，原始 Landsat 数据会进行标准化处理：

  • 热波段重建：在 Landsat 8 TIRS 热波段降级或未校准的运行期间，系统会合成占位符零方差热通道，以满足自动化 Landsat 云得分算法的内部接口要求。
  • Landsat 4/5 的合成全色生成：Landsat 4 和 5 TM 传感器缺少光学全色通道（波段 8）。为每个场景构建合成伪全色波段：
$$ I_{\text{raw45}}(Pan) = \text{mean}(I_{\text{raw45}}(RGB)) $$

下面 \(I_{\text{raw}}(\cdot)\) 是一个提供原始像素的图片函数，可能需要一个或多个波段作为参数。例如，\(I_{\text{raw}}(Pan)\) 返回原始全色波段，而\(I_{\text{raw}}(RGB)\) 返回原始 RGB 波段。mean() 是平均值函数。

校准为大气表观 (TOA) 反射率

使用 ee.Algorithms.Landsat.TOA 将 Landsat 原始场景处理为 TOA。

ee.Algorithms.Landsat.TOA 是一种图像函数，用于将 Landsat 原始数据转换为 Landsat TOA 数据 (\(I_{\text{TOA}}\))。

多传感器质量和云掩码

轨道幅宽边界渐变

距离衰减软 Alpha 遮罩有助于最大限度地减少 Landsat 轨道上的硬接缝：

$$ \text{Mask}_{\text{edge}} = \min\left(\text{Mask}_{\text{MODIS}}, \left(\text{Gaussian}_{6\text{km}}\left(\text{Mask}_{\text{raw}}\right)\right)^3\right) $$

其中， \(\text{Mask}_{\text{MODIS}}\) 是从 MODIS 地表反射率年度合成数据 (MOD09GA/MYD09GA) 派生的掩码， \(\text{Mask}_{\text{raw}}\)是 Landsat 原始场景的掩码，min() 是最小化函数， \(\text{Gaussian}\) 是具有指定标准差的高斯空间卷积。

Landsat 7 SLC 关闭水体掩码

2003 年 5 月，Landsat 7 上的扫描线校正器 (SLC) 发生故障，导致 Landsat 7 场景包含线性数据缺口，这可能会导致开阔水域出现严重伪影。使用模糊的 MODIS 水体掩码遮盖开阔水域上的 Landsat 7 像素，并根据镶嵌图制作年份优先使用 Landsat 4/5/8/9 或 Sentinel-2。

Landsat 云得分和异常恢复

使用 ee.Algorithms.Landsat.simpleCloudScore 为 Landsat 场景提供每像素 \(\text{Cloud}\) 得分（以 [0, 100] 为单位）。应用了两个额外的细化阶段：

  1. HSV 空间中未检测到的云恢复：通过将 \(I_{\text{TOA}}(RGB)\) 转换为色调-饱和度-值 (HSV) 空间并更新 \(\text{Cloud}\) 掩码，检测到接收到错误云得分（云 = 0）的高反射率、低饱和度云异常：
$$ \text{hsvCloud} = \begin{cases} 100, & \text{if } (\text{Cloud} = 0) \land (\text{S} \lt 0.3) \land (\text{V} \gt 0.5) \\ 0, & \text{otherwise} \end{cases} $$
$$ \text{Cloud} = \text{Cloud} + \text{hsvCloud} $$
  1. 自适应气候学阈值：用于舍弃多云像素（即 \(\text{Cloud} \gt \text{Threshold}_{\text{cloud}}\)）的阈值根据预先计算的空间云统计信息动态计算得出：
$$ \text{Threshold}_{\text{cloud}} = \text{Cloud}_{\text{p25}} + 10 $$

其中， \(\text{Cloud}_{\text{p25}}\) 是每个像素的第 25 百分位得分，该得分被限制在 25 到 65 之间。\(\text{Cloud}\) 对于经过验证的明亮陆地表面（沙漠沙丘、盐滩），阈值会提高：

$$ \text{Threshold}_{\text{cloud}} = \text{Threshold}_{\text{cloud}} + 10 \quad \text{if } (I_{\text{TOA}}(Pan) \gt 0.7) \land (\max(I_{\text{TOA}}(RGB)) \gt 0.6) $$

Sentinel-2 云得分

Sentinel-2 协调场景与 GOOGLE/CLOUD_SCORE_PLUS/V1/S2_HARMONIZED 配对：

  • 像素过滤：当晴空置信度 cs ≥ 0.60 且晴空累积概率 cs_cdf ≥ 0.70 时保留。

全色归一化和结构融合

为了利用 Landsat 7/8/9 的 15 米全色分辨率，同时避免传感器转换期间出现光谱失真，我们对 TOA 数据应用了全色归一化：

$$ \text{Gain}_{\text{3x3}} = \frac{\mu_{\text{3x3}}(\text{mean}(I_{\text{TOA}}(RGB))}{\mu_{\text{3x3}}(I_{\text{TOA}}(Pan))} $$
$$ \text{I}_{\text{sharpened}} = \frac{I_{\text{TOA}}(RGB)}{\text{mean}(I_{\text{TOA}}(RGB))} \times (I_{\text{TOA}}(Pan)) \times \text{Gain}_{\text{3x3}} $$

其中 \(\mu_{\text{3x3}}\) 是 3x3 像素邻域平均值函数。

针对 MODIS BRDF 基准进行辐射归一化

非线性对比度缩放：通过灰度系数调整将大气层顶部反射率缩放到 8 位显示范围 ([0,255])：

$$ I_{\text{scaled}} = \text{Visualize}\left(I_{\text{sharpened}}(RGB), \text{min}=0.02, \text{max}=0.50, \gamma=1.7\right) $$

其中， \(\text{Visualize}\) 是 ee.Image.visualize() 函数。除非另有说明，否则下文均指 8 位 RGB 颜色空间。

过度校正拒绝掩码：防止过度平滑，从而保护真实的快速地表覆盖变化：

$$ \text{Mask}_{\text{valid}} = (|I_{\text{scaled}}(R) - I_{\text{BRDF}}(R)| \le 60) \land (|I_{\text{scaled}}(G) - I_{\text{BRDF}}(G)| \le 40) $$

其中， \(I_{\text{BRDF}}\) 可根据 MODIS BRDF 调整后的反射率 (MCD43A4) 复合数据生成 8 位 RGB 数据。

低通空间场提取：

$$ \Delta_{\text{spatial}} = \text{Gaussian}_{20\text{km}}\left((I_{\text{scaled}} - I_{\text{BRDF}}) \times \text{Mask}_{\text{valid}}\right) $$

雾度减去和反射率协调：

$$ I_{\text{normalized}} = I_{\text{scaled}} - \Delta_{\text{spatial}} $$

陆地空间锐化

将拉普拉斯边缘卷积应用于 \(I_{\text{scaled}}\)和 \(I_{\text{normalized}}\) 的组合，仅在陆地表面上应用。

$$ I_{\text{sharp}} = \text{Sharpen}(I_{\text{normalized}}, I_{\text{scaled}}) $$

时间中值合成

对于每个日历年 \(Y\)，使用 16 位精度的中值归约器将归一化且经过云掩码处理的影像集合归约为单个多波段合成影像：

$$ I_{\text{annual}}(Y) = \operatorname{median}_{t \in Y}\left(I_{\text{sharp}}(t)\right) $$
$$ \text{Cloud}_{\text{annual}}(Y) = \operatorname{median}_{t \in Y}\left(\text{Cloud}(t)\right) $$

生成的图片表示原始的年度镶嵌图，包含 red、green、blue 和 cloud 波段，以及未插值的有效像素掩码。请注意，时间 (\(t\)) 是图像函数的形参。

后处理：全局间隙填充和水体遮罩资源

为了生成交互式 Earth 延时摄影查看器和年度拼接图中使用的无缝、全球完整的基本地图，需要对原始年度合成图像进行后处理。

时间线性回归缺口填充（1999 年之前的镶嵌）

在 1999 年 Landsat 7 发射之前，全球卫星采集受到历史下行链路基础设施、机载磁带记录器限制和持续云覆盖的制约。因此，1999 年之前的年度镶嵌影像包含大量空间空白，尤其是在中非、东南亚、西伯利亚和亚马逊地区。

为了消除令人分心的灰色空白，同时保留时间过渡，该流水线实现了逐像素的时间线性回归插值：

对于每个目标年份 Y，该算法会在数十年期集合中搜索之前 (\(t_{\text{before}} \le Y\)) 和之后 (\(t_{\text{after}} \ge Y\)) 最新的有效非掩码像素观测结果。

线性回归拟合

普通最小二乘回归模型根据每个像素进行评估，其中边界时间观测结果包含自变量 [1, t] 和因变量 [Red, Green, Blue]：

$$ I_{\text{sharp}}(RGB, t) = \mathbf{m} \cdot t + \mathbf{b} $$

模拟价值评估

在目标年份 Y 评估线性轨迹：

$$ I_{\text{interp}}(RGB, Y) = \mathbf{m} \cdot Y + \mathbf{b} $$

观察叠加

将 Y 年的年度像素与插值像素拼接在一起：

$$ I_{\text{gapfilled}} = \operatorname{Mosaic}\left(I_{\text{interp}}, I_{\text{annual}}\right) $$

这样可确保保留相应年份的所有真实卫星观测数据，同时对历史数据缺口进行时间插值。

高纬度冰川填充

Landsat 轨道倾角限制了在极地附近（高于 82.6°N）的观测。对于缺乏光学覆盖的格陵兰北部内陆地区和北极冰架，系统会混合使用经过校准的归一化多年基准，以保持清晰、无缝的极地基础地图。

全球海洋测深和水体掩膜

在原始的年度镶嵌图中，海洋表面包含阳光反射、云阴影和瞬时波伪影。最终资源将开阔的海洋水域替换为全球阴影浮雕地形地貌底图：

  1. 海深阴影：NOAA ETOPO1 模型中的地形海拔高度使用海洋深度颜色渐变进行样式设置：
    • 深海深渊平原（-5000 米）：深海军蓝 (#000927)
    • 大陆坡（-1000 米）：板岩蓝 (#000E3A)
    • 大陆架（-100 米）：天蓝色 (#000E3B)
    • 海岸线（0 米）：皇家钴蓝 (#001146)
  2. 山体阴影和高斯高通增强：结合了分析性山体阴影 (ee.Terrain.hillshade) 和 5,000 米高斯非锐化遮罩，以突出显示海沟、大洋中脊和海山。
  3. 陆地/水体边界协调：
    • 将 Hansen 全球森林变化数据 datamask（区分陆地和海洋）与 MODIS 水掩膜 MOD44W 相结合。
    • 绘制内陆湖泊和封闭海域（例如，里海、五大湖、贝加尔湖、咸海），以保留自然水色动态。
    • 应用 3 阶段形态圆形内核缩减和 3,000 米高斯模糊，以平滑地缩小浅海沿岸水域与近海测深之间的边界，而不会出现硬剪切。

全局色彩平衡和局部对比度增强 (LCE)

  1. HSV 值提升和伽玛平衡：通道伽玛（$\gamma_R = 0.98、\gamma_G = 1.00、\gamma_B = 1.04$），然后在 HSV 色彩空间中将值通道提升 10%。
  2. 极地冰白点校正：标记高反射率极地表面（格陵兰、南极洲、高山冰盖），并将其平衡为纯中性白色 ([255, 255, 255])。
  3. 多尺度 LCE：通过拉普拉斯卷积、平均值过滤和高斯锐化掩模的组合，增强不同地貌类型的对比度和清晰度。
功能 原始年度镶嵌 全球总决赛年度马赛克
主要应用场景 科学分析、来源跟踪、机器学习训练 可视化底图、视频导出、全球延时摄影探索
空间完整性 陆地（存在数据缺口） 100% 全球完整性
Pixel 丢失处理 已遮盖（透明无数据） 线性回归插值
1999 年之前的热带空洞 以无数据间隙的形式保留 在时间上平滑插值
海洋水表示法 已遮盖或未归一化的水 样式化 NOAA ETOPO1 水深测量
包含的光谱频段 红色、绿色、蓝色、云端 (QA) 红色、绿色、蓝色
位深 8 位无符号整数 (0–255) 8 位无符号整数 (0–255)
Pixel Resolution：1984 年至 2014 年 30.0 米/像素 (EPSG:3857) 30.0 米/像素 (EPSG:3857)
Pixel 分辨率：2015 年至今 每像素 19.11 米 (EPSG:3857) 每像素 19.11 米 (EPSG:3857)
网格维度（1984 年至 2014 年） 1,335,834 × 1,198,340 像素 1,335,834 × 1,198,340 像素
网格维度（2015 年至今） 2,097,152 × 1,881,297 像素 2,097,152 × 1,881,297 像素

表 2. 原始地球延时摄影集与全球最终地球延时摄影集之间的详细技术比较。

局限性和分析注意事项

  1. 双分辨率时代（30 米与 19.11 米）：进行数十年时间序列分析的用户必须考虑 2015 年的分辨率过渡。1984 年至 2014 年的合成图像以每像素 30.0 米的网格化分辨率呈现，而 2015 年以后的图像由于集成了 Sentinel-2 MSI 数据，因此以每像素 19.11 米的网格化分辨率呈现。
  2. 插值像素（1999 年之前）：在“全球最终版”集合中，1999 年之前的镶嵌影像中缺失的像素会进行时间插值。在多年间隔期间，如果某个区域正在经历突然的土地利用转变（例如，快速森林砍伐或水库建设），则插值像素将描绘出逐渐的线性转变，而不是突兀的离散事件。
  3. 更改波段比率：辐射归一化可根据 MODIS BRDF 目标值优化各个场景之间的视觉一致性。虽然相对空间模式得以保留，但派生的光谱指数（例如 NDVI、EVI）与 2 级地表反射率值不同。
  4. 物候混合：高纬度北部地区（>60°N）代表仲夏条件（DOY 150-270），而温带和热带地区代表年度综合中位数。

归因

知识共享许可
此数据集已获得知识共享署名 4.0 国际版许可，并且需要提供以下提供方信息： 
Google Earth Timelapse (Google, Landsat, Copernicus)
包含经过修改的 Copernicus Sentinel 数据 [2015 年至今]。请参阅 Sentinel 数据法律声明。