About Earth Timelapse

地球縮時攝影資料集提供 1984 年至 2022 年的地球變化視覺記錄，這項資料集是由 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 Mercator (EPSG:3857) 中格狀排列。
  • 2015 年至今 (19.11 公尺 / 像素)：透過融合 ESA Copernicus Sentinel-2A/2B 多光譜儀 (MSI) 觀測資料 (10 公尺 / 20 公尺原生) 與 Landsat 8 和 9 觀測資料合成。圖片會以 Web Mercator (EPSG:3857) 網格顯示，解析度約為 19.11 公尺。

年度影像集涵蓋全球所有陸地和沿海地區，緯度範圍約為南緯 82.6° 到北緯 83.69°。要綜合分析超過 40 年的地球觀測資料，必須協調不同衛星群在運作壽命、感應器模式、軌道幾何和氣氛擾動方面的差異。

Earth Timelapse 管道會從 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)，可提升視覺效果。

來源資料和感應器規格

地球縮時影片集錦採用下列太空地球觀測任務的資料：

衛星和感應器規格

衛星任務 感應器 作業日期範圍 空間解析度 (MS / 全色) 擷取的頻譜帶 Earth Engine 集合 ID 馬賽克解析度時代
Landsat 4 主題對應器 (TM) 1982-07-01 – 1993-12-14 30 公尺/不適用 藍色、綠色、紅色、近紅外線、短波紅外線 1、熱能、短波紅外線 2 LANDSAT/LT04/C02/T1 30 公尺 (1984 年至 1993 年)
Landsat 5 主題對應器 (TM) 1984 年 3 月 1 日至 2012 年 12 月 31 日 30 公尺/不適用 藍色、綠色、紅色、近紅外線、短波紅外線 1、熱能、短波紅外線 2 LANDSAT/LT05/C02/T1 30 公尺 (1984 年至 2012 年)
Landsat 7 增強型主題製圖儀 Plus (ETM+) 1999 年 4 月 15 日至 2013 年 8 月 31 日 30 公尺 / 15 公尺 藍色、綠色、紅色、近紅外線、短波紅外線 1、熱感應、短波紅外線 2、全色 LANDSAT/LE07/C02/T1
LANDSAT/LE07/C02/T2		 30 公尺 (1999 年至 2013 年)
Landsat 8 操作陸地成像儀 (OLI) / 熱紅外線感應器 (TIRS) 2013 年 2 月 11 日至今 30 公尺 / 15 公尺 沿海、藍色、綠色、紅色、近紅外線、短波紅外線 1、短波紅外線 2、全色、卷雲、熱 LANDSAT/LC08/C02/T1_RT_TOA
LANDSAT/LC08/C02/T2_TOA		 30 公尺 (2013 年至 2014 年)
19.11 公尺 (2015 年後)
Landsat 9 操作陸地成像儀 2 (OLI-2) / TIRS-2 2021 年 11 月 1 日至今 30 公尺 / 15 公尺 沿海、藍色、綠色、紅色、近紅外線、短波紅外線 1、短波紅外線 2、全色、卷雲、熱 LANDSAT/LC09/C02/T1_TOA 19.11 公尺 (2021 年以上)
Sentinel-2A / 2B 多光譜儀 (MSI) 2015 年 6 月 23 日至今 10 公尺 / 20 公尺 / 60 公尺 藍色、綠色、紅色、紅邊 (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. 漢森全球森林變化 (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° 到 83.69°)：
    • 依「一年中的第幾天」(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 大氣層頂端資料 (\(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-off 水域遮罩

2003 年 5 月，Landsat 7 的掃描線校正器 (SLC) 故障，導致 Landsat 7 場景含有線性資料空白，可能造成開放水域出現嚴重影像瑕疵。開放水域的 Landsat 7 像素會使用模糊的 MODIS 水遮罩遮蓋，並視影像拼接製作年份，優先使用 Landsat 4/5/8/9 或 Sentinel-2。

Landsat 雲端評分和異常復原

系統會使用 ee.Algorithms.Landsat.simpleCloudScore 為 Landsat 場景提供每像素 \(\text{Cloud}\) 分數 (以 [0, 100] 為單位)。系統會套用兩個額外的精修階段：

  1. HSV 空間中未偵測到的雲層復原：高反射率、低飽和度的雲層異常會收到錯誤的雲層分數 (雲層 = 0)，但透過轉換 \(I_{\text{TOA}}(RGB)\) 為色調-飽和度-值 (HSV) 空間並更新 \(\text{Cloud}\) 遮罩，即可偵測到這類異常：
$$ \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 個百分位數，\(\text{Cloud}\) 分數則介於 25 到 65 之間。針對已驗證的明亮地表 (沙漠沙丘、鹽灘)，門檻會提高：

$$ \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 配對：

  • Pixel 過濾：當晴空信心指數 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])，並調整 Gamma 值：

$$ 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\)) 是圖像函式的參數。

後續處理：全球填補缺口和遮罩水域的資產

為了製作互動式 Google 地球縮時攝影檢視畫面和年度影像集中使用的無縫全球完整基本地圖，我們對原始年度複合圖像進行後續處理。

時間線性迴歸間隙填補 (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 公尺)：Azure Navy (#000E3B)
    • 海岸線 (0 公尺)：皇家鈷藍色 (#001146)
  2. 陰影和高通高斯增強：結合分析陰影 (ee.Terrain.hillshade) 和 5000 公尺高斯銳利化遮罩，強調海溝、中洋脊和海山。
  3. 陸地/水域界線協調：
    • 結合 Hansen 全球森林變化 datamask (區分陸地和海洋) 與 MODIS 水體遮罩 (MOD44W)。
    • 繪製內陸湖泊和封閉海域 (例如：裡海、五大湖、貝加爾湖、鹹海)，以保留自然水色動態。
    • 套用 3 階段的形態圓形核心縮減和 3000 公尺的高斯模糊，平滑地縮減淺水海岸水域和近海深水地形之間的邊界，不會出現硬性剪裁。

全域色彩平衡和局部對比度增強 (LCE)

  1. HSV 值提升和 Gamma 平衡：通道 Gamma ($\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 Resolution：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 年之後的影像則以每像素 19.11 公尺的網格呈現，這是因為整合了 Sentinel-2 MSI 資料。
  2. 內插像素 (1999 年前)：在「全球最終」集合中，1999 年前鑲嵌圖像中缺少的像素會經過時間內插。在經歷多年間隔的期間，如果某個區域的土地使用情況突然發生變化 (例如快速砍伐森林或建造水庫)，內插像素會呈現漸進式線性變化，而非突然的離散事件。
  3. 改變波段比率：輻射歸一化會根據 MODIS BRDF 目標，針對各場景最佳化視覺一致性。雖然相對空間模式會保留，但衍生光譜指數 (例如 NDVI、EVI) 與第 2 級地表反射率值不同。
  4. 物候混合：高緯度北部地區 (60°N 以上) 代表仲夏條件 (DOY 150–270)，而溫帶和熱帶地區則代表年度綜合中位數。

歸因

創用 CC 授權
本資料集採用創用 CC 姓名標示 4.0 國際授權，並須註明出處： 
Google Earth Timelapse (Google, Landsat, Copernicus)
包含 [2015 年至今] 經修改的 Copernicus Sentinel 資料。請參閱「Sentinel 資料法律聲明」。