1985-2023年无云Landsat年度合成影像数据集：技术实现与应用指南

1次阅读

没有评论

共计 2660 个字符，预计需要花费 7 分钟才能阅读完成。

1985-2023 年无云 Landsat 年度合成影像数据集填补了长时序、全覆盖、高质量遥感数据的空白。其主要价值体现在：

长时间跨度：覆盖近 40 年地表变化，支持长期环境监测
无云干扰：通过合成算法消除云层影响，提升数据可用性
一致性处理：统一校正的反射率数据保证时序可比性

典型应用包括：

土地利用 / 覆盖变化检测（如森林砍伐监测）
农作物生长状况年度评估
城市扩张分析
生态环境演变研究

处理全球范围、多时相的 Landsat 数据（单景约 1GB）需要：

采用分块存储策略（如 COG 格式）
建立金字塔索引加速空间查询
使用云存储 + 本地缓存混合架构

并行计算框架选择（Dask/Spark）
内存优化技术（分块处理 / 数据压缩）
分布式任务调度（如 Kubernetes 集群）

严格的辐射定标与大气校正
多时相数据配准（亚像元级精度）
异常值检测与修复

数据获取：从 USGS EarthExplorer 或 Google Earth Engine 批量下载
质量筛选：基于 QA 波段剔除低质量像元
辐射归一化：应用 FLAASH 大气校正模型
几何校正：使用 DEM 数据消除地形畸变

采用改进的 FMask 算法，主要步骤：

# 基于 SCL 分类的云检测示例
import numpy as np

def detect_clouds(scl_band):
    """
    基于 Landsat SCL 分类结果检测云层
    :param scl_band: 场景分类波段(0-11)
    :return: 云掩膜(1= 云,0= 非云)
    """
    cloud_mask = np.isin(scl_band, [3,8,9,10])  # 3= 云影,8= 中云,9= 高云,10= 薄卷云
    cirrus_mask = scl_band == 10
    # 薄云增强处理
    return cloud_mask | (cirrus_mask & (np.random.random(scl_band.shape) > 0.7))

采用最佳像素合成 (BAP) 算法：

计算各像元在所有可用影像中的得分：
云量权重（0-1）
植被指数权重（NDVI）
传感器观测角度权重
选择综合得分最高的观测值
时空插值填补缺失数据

import rasterio
import xarray as xr

def load_annual_composite(year):
    """加载年度合成数据"""
    with rasterio.open(f'LC08_{year}_BAP.tif') as src:
        return xr.DataArray(src.read(),
            dims=('band','y','x'),
            coords={'band': ['B2','B3','B4','B5','B6','B7']}
        )

# 示例：加载 2015 年数据并计算 NDVI
ds_2015 = load_annual_composite(2015)
red = ds_2015.sel(band='B4')
nir = ds_2015.sel(band='B5')
ndvi_2015 = (nir - red) / (nir + red)

from sklearn.ensemble import RandomForestClassifier

def detect_landcover_change(year1, year2):
    """基于随机森林的变化检测"""
    # 1. 加载两年数据
    img1 = load_annual_composite(year1).stack(pixel=('y','x'))
    img2 = load_annual_composite(year2).stack(pixel=('y','x'))

    # 2. 计算差异特征
    diff_features = np.concatenate([(img2 - img1).values.T,  # 波段差值
        (img2 / (img1+1e-6)).values.T  # 波段比值
    ], axis=1)

    # 3. 训练分类器（示例简化，实际需标注数据）clf = RandomForestClassifier(n_estimators=100)
    # 假设已有训练数据 X_train, y_train
    clf.fit(X_train, y_train)

    # 4. 预测变化类型
    return clf.predict(diff_features).reshape((img1.y.size, img1.x.size))

分块处理：将大数据分块处理避免内存溢出

# Dask 分块示例
import dask.array as da
big_image = da.from_zarr('landsat_annual.zarr', chunks=(1, 2048, 2048))

延迟加载：只在需要时读取数据

with rasterio.open('large.tif') as src:
    window = Window(0, 0, 1024, 1024)
    subset = src.read(window=window)

并行计算：使用多进程处理不同区域

from concurrent.futures import ProcessPoolExecutor

def process_tile(tile):
    return heavy_computation(tile)

with ProcessPoolExecutor() as executor:
    results = list(executor.map(process_tile, tile_list))