简介面向遥感与地理信息领域的 Python 开发者资源提供了 FSDAF 时空融合算法的完整工程实现帮助学习者解决多时相、多源遥感影像在时间与空间分辨率上的综合增强问题适用于地表覆盖变化监测、作物长势分析、城市扩展研究等场景。压缩包共计 361 个文件、约 7.72MB其中包含 311 个 Python 源程序覆盖从数据读取、预处理到融合计算与结果导出的全流程另有少数配置文件、模型权重、批处理脚本与说明文档便于直接运行和二次开发。包内还随附 Python 虚拟环境相关文件可快速搭建依赖环境。目前已有 1140 人学习下载适合具备一定 Python 和遥感基础、希望深入理解 FSDAF 算法原理或将其迁移到自身研究数据的进阶使用者。通过研读源码与文档能够掌握时空融合的完整处理链条、调参思路与评价方法为后续科研或工程实践提供可直接参考的代码框架。1. FSDAF遥感影像时空融合为什么说它是“性价比最高”的融合方向在做地表覆盖变化、作物长势监测这类活儿时手里常常只有两种数据一种像 MODIS每天都能过境但一个像元几百米种了几种作物根本分不清另一种像 Landsat空间分辨率够用但好几天甚至十几天才能拿到一景遇到连续阴雨就彻底断档。FSDAF遥感影像时空融合要解决的就是这个矛盾——用算法把这两堆影像合成出一批“既高分辨率、又高时间密度”的新影像。FSDAF 的英文全称是 Flexible Spatiotemporal DAta Fusion这个名字里的 Flexible 很关键它对地表异质性、物候变化和突变的适应能力比早期方法强不少目前也是主流时空融合方案里最容易用 Python 落地的一个。这篇文章面向的是手里已经有配准好的 MODIS 和 Landsat 影像想把融合代码跑起来、并且能看懂每个环节在干什么的从业者。我会从 FSDAF “预测 残差修正”的底层逻辑说起然后给出一套完整的 Python 实现结构与核心模块代码最后把我实际跑数据时踩过的坑和验证方法一并写出来。不保证代码贴上去就能出论文级结果但能保证每条参数你都调得有依据。2. FSDAF 是怎么把两堆影像“算”成第三堆的核心机制与计算代价2.1 时空融合的输入输出关系先搞清楚你要准备几组影像所有时空融合方法都建立在同一个假设上粗分辨率影像比如 MODIS和细分辨率影像比如 Landsat在同一个时刻观测的是同一片地表它们之间的差异只来自空间分辨率而不是来自地物发生了变化。这句话看起来简单但决定了整个算法的输入输出结构。FSDAF 需要至少三组数据一组是预测时刻 t2 的粗分辨率影像一组是已知时刻 t1 的粗分辨率影像还有一组是 t1 时刻的细分辨率影像。输出则是 t2 时刻的细分辨率影像。如果你手里有 t1 时刻的细分辨率真实影像还可以把 t2 时刻的高分辨率预测结果和真实影像做精度验证——这是判断融合效果最直接的办法。为了让它跑得稳输入影像不能是随便裁剪的。我在第一次跑的时候用的是覆盖某农田试验区的一景 Landsat-8 和一景 MODISLandsat 经过辐射定标和大气校正MODIS 用的也是地表反射率产品。两组数据必须做三件事一是投影坐标系完全一致二是像元边界严格对齐三是数据范围大致覆盖同一块区域。MODIS 像元是 500 米Landsat 是 30 米比例大约是 16.67 倍所以一个 MODIS 像元大概对应 17×17 个 Landsat 像元这个比例在后面的窗口参数设置里会反复用到。2.2 两次预测与残差分配FSDAF 内部到底走了哪几步FSDAF 的核心可以拆成“两次预测 一次残差修正”这也是它区别于早期 STARFM时空自适应反射率融合模型的地方。STARFM 的做法是把邻域内相似像元的权重算出来直接加权而 FSDAF 则先估计地表覆盖类型的变化再做时序外推最后用残差把细节补回去。第一步算法对 t1 时刻的细分辨率影像做分类或聚类把像元划分成若干类别比如水体、裸土、作物、林地。然后在每个类别内部建立 t1 和 t2 两个时相粗分辨率影像之间的线性回归关系得到每个类别的变化趋势。为什么要在类别内做回归因为不同地类的物候变化节奏完全不同混在一起回归会把变化趋势平均掉比如水体几乎不变作物正在快速生长混在一起的结果就是水体被高估了变化、作物被低估了变化。第二步算法用 t1 时刻的细分辨率影像结合回归预测的结果生成一个 t2 时刻的初步融合结果。但这一步的误差很大因为它假设每个细像元的变化规律只取决于它属于哪个类别而忽略了类别内部的空间异质性。所以就有了第三步计算 t1 和 t2 时刻粗分辨率影像之间的真实差值减去预测得到的差值剩下的部分就是残差。这个残差被空间插值分配到每个细像元上再做一次修正得到最终的融合结果。用一句话概括就是先用类别回归搭骨架再用粗分辨率的变化残差补血肉。2.3 为什么 FSDAF 不挑地类关于“灵活性”的一个解释早期方法对环境变化特别敏感。比如某块地在 t1 和 t2 之间发生了突变——洪水淹没、森林火灾、城市扩张STARFM 这类基于局部加权的方法往往会把突变“抹匀”因为它的权重机制倾向于依赖邻近相似像元而突变区域的邻近像元还是老样子。FSDAF 之所以在这种场景下更稳是因为它的残差修正机制没有依赖“邻近像元相似”这个强假设而是直接计算全局粗分辨率影像的差值再把这些差值插值回细分辨率网格。这意味着即使某个区域的地表类型发生了剧烈变化只要粗分辨率影像捕捉到了这个变化比如 MODIS 影像上能看到洪水前后反射率的明显差异FSDAF 就能通过残差分配把这个信号带回细分辨率结果里。代价是计算量明显增大分类、回归、插值、迭代每一步都在全图上做矩阵运算而且插值环节用的是薄板样条数据量大时内存会吃紧。我个人的习惯是如果研究区地形破碎、地类复杂FSDAF 是首选如果是大范围均质的农田或草地STARFM 也能凑合但 FSDAF 的结果通常更细腻尤其是在边界区域。3. 把 FSDAF 写成 Python工程架构与最小复现目录3.1 一个顺手的文件组织别把代码全塞进一个脚本里网上能搜到的 FSDAF Python 实现大多散落在各个论文复现项目或个人维护的仓库里代码风格差异很大有的把所有逻辑塞进一个脚本有的拆成十几个模块却没有一处注释。我的建议是自己重写一套结构清晰、参数透明出了问题也知道去哪找。我习惯的文件组织如下fsdaf_project/ ├── data/ │ ├── landsat_t1.tif │ ├── modis_t1.tif │ └── modis_t2.tif ├── config.yaml ├── utils.py ├── fsdaf_core.py ├── main.py └── output/这里 config.yaml 存所有可调参数utils.py 负责数据读写和预处理fsdaf_core.py 是算法主体main.py 是入口脚本。这样拆分的好处是算法逻辑不受 IO 干扰参数调整不需要改代码换数据集时只动配置文件和 data 目录。3.2 数据结构与参数管理YAML 配置比硬编码好一百倍FSDAF 的参数不算多但每个都对结果有直接影响。我把它们统一放到 config.yaml 里每次跑实验前先看一遍配置避免在代码里翻找魔法数字。# config.yaml data: coarse_dir: ./data/modis fine_dir: ./data/landsat t1_date: 2023-05-01 t2_date: 2023-05-17 params: coarse_res: 500 # 粗分辨率影像像元大小米 fine_res: 30 # 细分辨率影像像元大小米 window_size: 1500 # 搜索窗口边长米通常取粗像元的整数倍 classes: 5 # 分类类别数不宜太少也不宜太多 max_iter: 3 # 迭代修正次数一般 1~3 足够 tps_smooth: 0.1 # 薄板样条插值的平滑参数参数说明window_size 建议不小于粗像元边长的 3 倍比如 MODIS 500 米分辨率窗口取 1500 米刚好覆盖 3×3 个 MODIS 像元太小会丢失局部变化信息太大则计算量爆炸。classes 取决于研究区地类复杂度5 到 10 是个合理区间超过 15 会导致某些类别像元数过少回归结果不稳定。max_iter 不是越大越好迭代到第三次以后残差基本收敛反而可能引入噪音。3.3 数据读取与预处理GDAL 是绕不开的依赖遥感影像读写绕不开 GDALPython 里最顺手的封装是 rasterio底层就是 GDAL。读取影像后要马上确认三个信息投影、像元大小、有效值范围。不一致的第一时间做重投影和重采样不要拖到后面才发现问题。# utils.py import numpy as np import rasterio from rasterio.warp import reproject, Resampling def read_geotiff(path): with rasterio.open(path) as src: data src.read(1).astype(np.float32) profile src.profile transform src.transform crs src.crs return data, profile, transform, crs def align_to_fine(coarse_data, fine_profile, coarse_transform, coarse_crs): 将粗分辨率影像重采样到细分辨率影像的网格。 dst_transform fine_profile[transform] dst_crs fine_profile[crs] dst_shape (fine_profile[height], fine_profile[width]) aligned np.zeros(dst_shape, dtypenp.float32) reproject( sourcecoarse_data, destinationaligned, src_transformcoarse_transform, src_crscoarse_crs, dst_transformdst_transform, dst_crsdst_crs, resamplingResampling.bilinear, ) return aligned逻辑说明read_geotiff 函数读入单波段影像并转成 float32因为后续计算涉及大量浮点运算int16 容易溢出或损失精度。align_to_fine 把 MODIS 重采样到 Landsat 的网格上注意这里用的是双线性重采样而不是最近邻因为后续的像元级别计算需要平滑过渡最近邻会产生块状效应。如果研究区范围大建议先读影像元数据确认行列数差异50 万行以上的数据要考虑分块处理不然 align 的时候内存会爆。3.4 主干流程脚本入口越简单越好main.py 只做一件事按顺序调用预处理、算法、输出中间打印每步耗时这样跑长任务时至少知道卡在哪了。# main.py import yaml import time import numpy as np from utils import read_geotiff, align_to_fine from fsdaf_core import fsdaf_predict if __name__ __main__: with open(config.yaml, r) as f: cfg yaml.safe_load(f) t0 time.time() fine_t1, fine_profile, fine_transform, fine_crs read_geotiff( f{cfg[data][fine_dir]}/landsat_t1.tif ) coarse_t1, coarse_profile, coarse_transform, coarse_crs read_geotiff( f{cfg[data][coarse_dir]}/modis_t1.tif ) coarse_t2, _, _, _ read_geotiff( f{cfg[data][coarse_dir]}/modis_t2.tif ) print(f数据读取完成耗时 {time.time() - t0:.1f}s) coarse_t1_aligned align_to_fine( coarse_t1, fine_profile, coarse_transform, coarse_crs ) coarse_t2_aligned align_to_fine( coarse_t2, fine_profile, coarse_transform, coarse_crs ) print(影像对齐完成) result fsdaf_predict( fine_t1fine_t1, coarse_t1coarse_t1_aligned, coarse_t2coarse_t2_aligned, cfgcfg[params], ) print(f融合计算完成总耗时 {time.time() - t0:.1f}s) with rasterio.open( output/fsdaf_fusion.tif, w, driverGTiff, heightresult.shape[0], widthresult.shape[1], count1, dtypefloat32, crsfine_crs, transformfine_transform, ) as dst: dst.write(result, 1)逻辑说明入口脚本中所有数据处理都走 utils 和 fsdaf_core不在入口里写任何算法逻辑。配置文件的读入放在最前面后面所有函数的参数都从 cfg 里取。输出文件直接写入 output 目录投影和仿射变换参数沿用 Landsat t1 的保证地理参考不乱。4. 从论文到代码核心计算模块的逐步实现4.1 高分辨率端的像元定位与分类映射先给地表分门别类FSDAF 的第一步是在细分辨率影像上做分类。常见实现里为了省事直接用 KMeans 聚类而不是用监督分类因为不需要训练样本而且聚类结果对回归来说是够用的。KMeans 的类别数由 config 里的 classes 控制。# fsdaf_core.py from sklearn.cluster import KMeans def classify_fine_image(fine_t1, n_classes): 对细分辨率影像做无监督分类返回每个像元的类别标签。 rows, cols fine_t1.shape # 展平成向量忽略 nodata 像元 valid_mask fine_t1 0 values fine_t1[valid_mask].reshape(-1, 1) kmeans KMeans(n_clustersn_classes, random_state42, n_init10) labels_flat kmeans.fit_predict(values) labels np.zeros(rows * cols, dtypenp.int16) labels[valid_mask.ravel()] labels_flat labels labels.reshape(rows, cols) # 类别编号从 1 开始0 留给无效像元 labels labels 1 labels[~valid_mask] 0 return labels参数说明n_init10 是必要的KMeans 对初始中心敏感n_init10 意味着跑 10 次取最好结果代价是时间变长但稳定性收益值得。valid_mask 用大于 0 来过滤是因为遥感影像的 nodata 值通常是 0 或负数如果数据里有水体且 DN 值很低注意确认它能不能被保留。分类做完后需要对每个类别统计细分辨率影像的均值以及在对应位置的粗分辨率影像均值。这为下一步回归准备了配对数据。def class_means(fine_t1, coarse_t1, labels, n_classes): 每个类别内细影像和粗影像的平均值。 class_fine_means np.zeros(n_classes 1) class_coarse_means np.zeros(n_classes 1) for c in range(1, n_classes 1): mask labels c if mask.sum() 0: class_fine_means[c] fine_t1[mask].mean() class_coarse_means[c] coarse_t1[mask].mean() return class_fine_means, class_coarse_means逻辑说明类别均值的作用是后续的推断基准。某个类别如果像元太少均值会很不稳定所以前面建议 classes 不要太大就是这个原因。4.2 粗分辨率端的回归预测变化趋势的数学表达这一步要做的是找到从 t1 到 t2 每个地类的变化规律。实现思路是利用 t1 和 t2 的粗分辨率影像计算每个类别对应的变化量然后用这个变化量去驱动细分辨率的预测。def coarse_change_per_class(coarse_t1, coarse_t2, fine_t1, labels, n_classes, scale500 // 30): 基于粗分辨率影像计算每个类别在 t1 和 t2 之间的反射率变化。 参数 scale 用于将粗像元划分为细像元块默认 500m/30m 取整。 实际实现中 coarse_t1 和 coarse_t2 已经重采样到细网格 因此直接用类别掩膜求均值即可。 class_delta np.zeros(n_classes 1) for c in range(1, n_classes 1): mask labels c if mask.sum() 0: delta coarse_t2[mask].mean() - coarse_t1[mask].mean() class_delta[c] delta return class_delta严格来说论文里这一步是用每个类别内部细像元上的粗分辨率变化做回归确定一个线性系数。工程实现上可以简化成“类别内均值差”因为多数情况下类别内反射率变化近似线性均值差已经能代表整体趋势。如果你的研究区物候差异特别大比如同一类别里有早熟和晚熟品种建议改成在每个类别内再对粗分辨率差值做一次局部线性回归斜率和截距分别存下来预测时逐个像元外推。4.3 残差分配薄板样条插值把细节补回来残差计算是 FSDAF 的精髓。用粗影像 t2 减去 t1 得到真实变化减去预测变化得到一个残差场。这个残差场的空间分辨率是粗分辨率级别的需要插值到细分辨率。工程上常用薄板样条thin plate spline或者普通克里金但考虑到纯 Python 实现的性能我更推荐用 scipy 的 griddata 配合 cubic 插值作为近似替代。我个人测试下来薄板样条的平滑性更好残差场的边缘过渡自然cubic 方法在像元数量大时两者差异不大但 cubic 快得多。如果追求和论文保持一致的算法对齐用 TPS如果在工程效率上有要求用 cubic。from scipy.interpolate import griddata def interpolate_residual(coarse_residual, fine_shape, coarse_grid): 将粗分辨率残差场插值到细分辨率网格。 rows, cols fine_shape fine_y, fine_x np.mgrid[0:rows, 0:cols] coarse_y, coarse_x coarse_grid # 粗网格的点位坐标 values coarse_residual.ravel() # 去掉无效值 mask np.isfinite(values) interpolated griddata( (coarse_x[mask], coarse_y[mask]), values[mask], (fine_x, fine_y), methodcubic, ) return interpolated参数说明coarse_grid 是粗分辨率影像每个像元的行列坐标注意这里用的是阵列坐标而不是地理坐标这样插值结果直接落到细网格的行列空间省去投影变换。cubic 插值在边界区域容易产生 overshoot也就是出现过高或过低的异常值如果结果里有明显跳变考虑用 methodlinear 或者对插值后的结果做一次截止值裁剪。4.4 迭代式修正与锐化输出把多个结果叠起来FSDAF 论文里的最终版本是做了两次预测第一次用类别回归给出基础预测第二次用残差修正。如果把这两个步骤叠起来跑多轮会得到更平滑的结果。工程里我会控制迭代次数在 2 到 3 次每轮用上一次的输出作为下一轮的 fine_t1 输入重新分类、回归、残差修正。def fsdaf_predict(fine_t1, coarse_t1, coarse_t2, cfg): FSDAF 主流程先做过分类回归再做残差修正最后迭代输出。 n_classes cfg[classes] max_iter cfg[max_iter] current_fine fine_t1.copy() current_coarse_t1 coarse_t1.copy() for it in range(max_iter): # 1. 分类 labels classify_fine_image(current_fine, n_classes) # 2. 每个类别的粗分辨率变化 class_delta coarse_change_per_class( coarse_t1, coarse_t2, current_fine, labels, n_classes, ) # 3. 基础预测细影像 t1 类别变化 base_pred current_fine.copy() for c in range(1, n_classes 1): mask labels c base_pred[mask] class_delta[c] # 4. 残差真实变化减去基于类别的预测变化 predicted_coarse_change coarse_change_per_class( coarse_t1, base_pred, labels, n_classes ) # 粗分辨率尺度上的残差 residual (coarse_t2 - coarse_t1) - predicted_coarse_change[labels] # 插值到细分辨率 residual_fine interpolate_residual( residual, fine_t1.shape, coarse_grid ) # 5. 修正 current_fine base_pred residual_fine return current_fine逻辑说明这里的代码不是逐行对照论文公式而是把论文的思路提炼成了可运行的工程版本。关键区别是论文在对每个粗像元做残差时考虑了邻域权重我这里直接把残差场做了插值效率高很多但边界细节会稍弱。如果你的计算资源充裕可以在 interpolation 前对残差场做一次高斯平滑消除条带效应再进入插值环节。参数说明max_iter 别设过大超过 3 次后会开始放大插值带来的噪音。cls n_classes 过低时类别变化估计会偏粗边缘会糊过高时小类别回归不稳定可能出现局部亮斑。我一般先从 5 类和 2 次迭代起步看一眼结果再调整方向。5. FSDAF Python 实现避坑指南从“能跑”到“结果能看”5.1 影像配准误差直接被放大现象融合结果出现明显的“重影”地物边缘有双重轮廓像是两张影像没叠准。原因MODIS 和 Landsat 不是同一颗卫星过境时间、传感器视角都不同如果只靠原始几何信息直接叠加偏移量可能达到 1 到 2 个 MODIS 像元。FSDAF 残差修正对这种偏移非常敏感因为残差是逐像元计算的配准误差会直接融入残差场。解决在预处理阶段必须做一次影像配准。常见做法是拿 Landsat 作为基准用 MODIS 的几何文件做二次配准或者在 Python 里利用影像上明显的道路交点、水体边界手动选取控制点。如果项目里允许推荐用高分辨率影像先做一次全局配准再裁剪研究区。5.2 反射率放大倍数没换算对现象融合结果整体偏亮或偏暗直方图和真实影像差一个倍数关系。原因Landsat 和 MODIS 虽然都叫“地表反射率”但产品缩放尺度不一样。有的数据给的是 0 到 1 的反射率有的给的是 0 到 10000 的整数FSDAF 的回归和差值运算里没统一量纲结果自然不对。解决在读数据之后和存结果之前分别做一次归一化检查。我的习惯是把所有影像统一到 0 到 1 的浮点反射率乘以系数后再进入算法输出时再乘回去。这里不能在预处理时省略因为差值和回归对常数偏移免疫但对缩放倍数非常敏感。5.3 水体云影区域产生极端异常值现象融合结果在水体边缘或云影覆盖区出现负值或远大于 1 的反射率。原因这些区域的粗分辨率影像混合了水体和陆地信号类别的均值差不能代表细像元的变化插值残差在这些边界处会产生 overshoot。解决对输出结果做物理范围裁剪反射率小于 0 的按 0 处理大于 1 的按 1 处理。更深一层是输入时就做云掩膜把这些区域标记成无效在分类和回归阶段就把它们排除在外。别指望算法自己处理异常像元它只会把异常扩散到邻域。5.4 窗口大小缩短导致结果出现“马赛克”现象融合结果在类别边界处出现明显的方块状痕迹如同打了马赛克。原因窗口参数取太小残差插值只参考了单个粗像元局部细节无法从邻域获得信息。FSDAF 在计算类别均值时如果只统计窗口内的像元窗口太小会让类别均值的方差变大。解决把 window_size 调大到粗像元边长的 3 倍以上同时检查粗分辨率影像是否被重采样到过高分辨率比如 500 米数据重采样到 10 米像元数量膨胀但信息量并没有增加插值就会在局部制造虚假细节。5.5 迭代次数增多后结果反而变差现象max_iter 从 2 调到 5融合结果不仅没有更精细反而出现周期性条纹。原因迭代过程里每一轮都会产生新的插值残差而这些残差中有一部分来自插值的人工痕迹不会随着迭代消失只会被逐轮放大。解决固定 max_iter2 或 3并用真实 t2 细分辨率影像做精度对比。如果两轮迭代后 RMSE均方根误差不再下降就不要再加了。可以把迭代次数写进配置跑一组对比实验用实际数据决定而不是靠感觉。6. 验证方法写完后再谈两个进阶技巧6.1 用时间序列验证让融合结果从“算出来”变成“能用”迭代次数和类别数确定之后我习惯做一组时间序列验证选取 t1 和 t2 之外的第三天 t3用 t1 和 t3 的粗分辨率影像再加 t1 细影像做融合得到 t3 的预测细影像与真实的 t3 细影像对比。对比指标用 RMSE 和结构相似度SSIM。RMSE 反映数值准确度SSIM 反映纹理结构保真度。如果 SSIM 低于 0.7说明空间细节损失偏大优先检查残差插值环节如果 RMSE 偏高优先检查类别数和回归稳定性。这个方法比只看一景图主观判断靠谱得多。6.2 从双时相到连续时间序列的扩展FSDAF 本身是双时相融合器但遥感应用里经常需要连续时间序列。我的做法是滑动窗口式调用比如有 1 月到 10 月每隔 8 天一景的 MODIS 和每个月一景的 Landsat就用每一景 Landsat 作为 t1 基准分别融合它到下一个月 Landsat 日期之间所有的 MODIS 日期。这样每段融合都只跨越不超过 30 天物候变化小融合误差可控。这个思路比自己硬套长时序版本靠谱得多。代码到这里能跑通、能调参、能验证就算真正把 FSDAF 用起来了。回头看这个项目我最大的教训是千万别上来就调算法细节先拿一组干净数据把流程跑通再回头看参数。概率再高的算法也救不了脏数据。希望这篇整理对你有点用少走两步弯路。本文还有配套的精品资源点击获取