简介面向遥感图像处理与时空融合研究者的Python实现资源围绕FSDAF算法提供从数据预处理到融合结果评估的完整代码流程适合地理信息科学、环境监测等领域具备一定Python基础的学习者参考。压缩包共361个文件主体为311个py脚本另含exe可执行程序、hdr头文件、xml配置、txt说明、pth模型权重、yaml参数等分别承担运行环境、参数配置、数据读取与算法模块等功能整体包体约7.72MB。目前已有1138人学习下载。借助Landsat与MODIS多时相示例数据可复现FSDAF融合实验并对照分类结果掌握融合策略的参数配置与实现细节代码模块划分清晰便于二次开发可用于土地覆盖变化分析、植被监测等应用场景。 做遥感时序分析的人大概都经历过这种纠结手里MODIS影像一个月能攒十来景时间序列画出来几乎连续可500米分辨率上连城市主干道和农田边界都分不清Landsat倒是能把地块看得清清楚楚但16天重访一次再碰上云层遮挡一个生长季能用上的无云影像可能就三到四景。FSDAFFlexible Spatiotemporal Data Fusion就是为了打破这个“时间-空间跷跷板”被提出来的。简单说它用MODIS的时间变化信息去补Landsat的稀疏观测最终输出一景30米分辨率、时间频率却跟MODIS看齐的融合影像。这篇文章不是论文翻译是我从准备数据到写完第一版FSDAF遥感影像时空融合python代码再到把结果定量评估完的完整记录。里面会讲清楚算法每一步在做什么、为什么这样做以及文档里找不到的那些坑。适合正在做植被物候、农田监测、城市扩张或灾害应急又不想被商业软件绑住手脚的人。1. 为什么需要FSDAFMODIS与Landsat的“时空跷跷板”困境1.1 两类数据的互补特性先看一组实际参数对比传感器重访周期空间分辨率优势劣势MODIS1-2天250m-1km时间密度极高云层筛完后仍有大量可用影像空间分辨率不足异质区域无法使用Landsat16天30m空间细节清晰适合地块级别分析时间间隔长有效无云观测太少Sentinel-25天10-20m空间时间兼顾2015年后才有数据无历史长序列时间重访和空间分辨率在物理上互相制约传感器设计很难两头占优。时空融合的思路就是用同一个时间段里大量存在的MODIS影像作为“高频信息”用两景Landsat影像分别锁定“起始状态”和“参考状态”把tk时刻的MODIS变化“注入”到细分辨率的Landsat图像上从而得到一个理论上同时具备高时间分辨率和高空间分辨率的tk时刻预测图。这项技术最典型的应用场景有三个物候过渡期内的作物生长监测、云干扰严重区域的地表反射率重建、以及历史数据缺失区域的Landsat序列补全。1.2 STARFM的局限与FSDAF改进点提到时空融合绕不开STARFM。STARFM假设从t1到t2之间地物端元光谱不变只通过寻找相似像元做加权预测。这个假设在均质农田区域效果不错但到了城市、山地、破碎耕地这类混合像元密集区STARFM很容易把一个像元里几种地物的变化混在一起导致地物边界模糊、局部细节丢失。FSDAF的核心改动是引入光谱线性解混和残差空间分配两个步骤。光谱线性解混把MODIS粗像元看成内部多种地物端元的线性混合先解出每一种土地覆盖类别在t1到tk之间的反射率变化量再把这个变化量按类别重新赋给细影像上的每一个像元。这样即便一个MODIS像元里面混杂了农田和裸土两者变化方向不同也能分别处理而不是像STARFM那样全部算平均。另一个改进是残差分配。即使光谱解混做得再准细预测结果上采样回粗分辨率后和真实MODIS观测之间仍然存在误差。FSDAF用薄板样条插值把这种粗尺度残差平滑地分配到每个细像元上补偿局部异质性带来的预测偏差。2. FSDAF核心流程拆解分类、端元、解混、残差分配2.1 算法五步走的整体架构FSDAF的通用输入是三景影像加上一景待预测影像信息t1时刻细分辨率影像如Landsat、t2时刻细分辨率影像用于验证做单时相预测时可不强制使用、t1时刻粗分辨率影像如MODIS、tk时刻粗分辨率影像待预测时刻。实际预测时最多用到t1和tk两个时刻的MODISt2主要是用来做精度验证。完整的预测流程可以归纳成五个步骤对t1细影像做非监督分类得到K个土地覆盖类别并计算每一类的端元光谱该类所有像元的平均反射率。基于t1细影像的类别比例对粗影像做光谱解混根据MODIS t1和tk影像的差异求解每一类的反射率变化量ΔF(k)。将类别级变化映射回细影像用t1细影像像元所属类别对应的ΔF(k)生成初始预测。计算粗尺度残差把初始预测结果上采样到粗分辨率网格和真实tk时刻MODIS相减得到残差面。用薄板样条TPS插值把残差分配到细分辨率网格叠加到初始预测上得到最终融合影像。很多开源版本还会在最后加一步“时相一致性修正”用t2时刻的Landsat做反向预测对上一步结果做加权平滑但这部分不是FSDAF的必选步骤第一次实现可以暂时不碰。2.2 关键公式与物理含义光谱解混是FSDAF的灵魂。一个MODIS粗像元在t1时刻的观测值可以近似表达为y(t1, i) Σ(k1..K) A(i,k) * E(t1, k) ε其中y(t1, i)是第i个粗像元的反射率A(i,k)是第k类在该像元内所占的面积比例E(t1,k)是第k类在t1时刻的端元光谱ε是噪声。A矩阵可以从t1细影像的分类结果按粗像元范围统计得到。如果假设端元光谱在t1到tk之间自身有变化而类别比例不变那么两边同时相减Δy(i) y(tk, i) - y(t1, i) Σ(k1..K) A(i,k) * ΔE(k)这里ΔE(k)就是第k类从t1到tk的反射率变化量也就是我后面代码里的change_by_class。实际求解时把每个波段分别处理用最小二乘法解一个带约束的线性方程组。残差分配用薄板样条插值数学上等价于找一个弯曲能量最小的光滑曲面让它尽量穿过所有粗网格控制点。比起普通的双线性插值薄板样条不会在像元边界上产生折痕能生成物理上更合理的平滑残差面这对异质区域尤其重要。3. 环境准备与数据预处理融合结果好不好七成看这里3.1 Python环境与依赖库FSDAF本身没有特别复杂的高性能计算需求纯Python加numpy就能实现实际跑起来的主要瓶颈是数据体的I/O和薄板样条插值的计算耗时。推荐直接用conda建一个干净的虚拟环境conda create -n fsdaf python3.10 conda activate fsdaf pip install numpy scipy scikit-learn rasterio GDAL四个核心库的分工numpy所有数组运算光谱解混的核心工具scikit-learn提供KMeans聚类用来给t1细影像做非监督分类scipy提供RBFInterpolator实现薄板样条插值rasterio读写GeoTIFF处理地理坐标与投影信息我不建议用osgeo的gdal直接做数组操作rasterio的接口更现代尤其在处理nodata和transform时体验好很多。3.2 Landsat/MODIS数据准备与对齐预处理阶段最重要的原则是一切以细影像网格为基准。我见过太多人直接在融合环节报尺寸不匹配一问才发现两套影像的角点坐标差了半个像元。具体分四步统一坐标系和网格把MODIS重投影到Landsat所在的UTM投影并用最近邻或双线性重采样到30米分辨率。MODIS原始反射率像元大小是500米MOD09GA或250米MOD09GQ重采样时务必把目标网格的角点和Landsat对齐。统一数值范围Landsat Collection 2 Surface Reflectance产品的DN值乘以0.0001才是反射率MODIS MOD09GA反射率同样需要乘以0.0001。两个产品原始量纲不一致不换算直接进算法结果会完全跑偏。云掩膜处理FSDAF对云和云阴影非常敏感。云区的反射率异常高会被算法当成真实光谱变化传递到融合结果里。建议根据Landsat的QA_PIXEL波段和MODIS的State_1km波段分别做掩膜对掩膜区域用同期的邻域像元做简单填补或在分类前直接排除。裁剪同一研究区用rasterio读入后再统一做window裁剪保证影像行列数一致。预处理做完建议先写一小段代码验证两幅影像的transform是否完全一致。一致的标准是分辨率、原点坐标、旋转参数一般都为0全部相等。这一步能省掉后续大量莫名其妙的bug。4. 核心Python代码逐段拆解分类、端元与残差分配的实现细节4.1 参数与数据装载我习惯把所有可调参数集中在一个地方方便反复测试。下面是数据装载的骨架代码import numpy as np import rasterio from sklearn.cluster import KMeans from scipy.interpolate import RBFInterpolator # 可调参数 N_CLASS 8 # 土地覆盖分类数 WINDOW_SIZE 9 # 相似像元窗口odd SCALE 50 # MODIS像元对应Landsat像元数500m/30m约等于17 # 实际用重采样后的行数比值 BAND_NUM 4 # 波段数 fine_t1_path landsat_t1.tif fine_t2_path landsat_t2.tif coarse_t1_path modis_t1.tif coarse_tk_path modis_tk.tif out_path fusion_tk.tif def read_tif(path): with rasterio.open(path) as src: return src.read(), src.transform, src.crs, src.nodata注意SCALE这个参数。如果你已经把MODIS重采样到了30米分辨率那SCALE理论上应该是1两幅影像尺寸完全一致。但这样做会引入重采样误差而且计算量成倍增加通常不建议。比较合理的做法是保留MODIS原始分辨率如500米在特征提取时按比例映射到Landsat网格上。实际代码里SCALE一般取行列数的比值写成变量方便后续处理。4.2 土地覆盖分类与端元光谱计算FSDAF不要求预先知道地物类别直接对t1细影像做聚类即可。我实测下来KMeans足够稳定ISODATA在这个场景下优势不明显但KMeans的n_init要设大一点避免局部最优。def get_class_and_endmember(fine_img, n_classes): h, w, bands fine_img.shape data fine_img.reshape(-1, bands) # 把nodata像元剔除避免污染聚类中心 valid np.all(np.isfinite(data), axis1) kmeans KMeans(n_clustersn_classes, random_state42, n_init10) labels_flat np.full(data.shape[0], -1, dtypeint) labels_flat[valid] kmeans.fit_predict(data[valid]) class_map labels_flat.reshape(h, w) endmember np.zeros((n_classes, bands)) for c in range(n_classes): mask class_map c if mask.sum() 0: endmember[c] fine_img[mask].mean(axis0) return class_map, endmember这里有个容易被忽略的细节必须剔除nodata像元再聚类。如果研究区里有水域或云层残留这些像元的光谱容易出现极端值拉偏聚类中心导致端元光谱失真。剔除后再聚类最后用np.full填充回原图保证空间位置不塌缩。端元光谱本质上是该类在某波段上的平均反射率FSDAF原论文称之为“初始端元”。后续光谱解混时这个端元会用来构建线性方程组所以它的准确性直接决定解混结果的质量。4.3 粗尺度光谱解混求解每类变化量这一步把MODIS t1和tk的差异分解成K个类别各自的变化量。对每个粗像元需要先统计其覆盖范围内的Landsat像元属于每一类的数量比例形成类别比例矩阵A然后用最小二乘解方程。def unmix_change(fine_t1, coarse_t1, coarse_tk, class_map, endmember, scale): h, w fine_t1.shape[:2] bands fine_t1.shape[2] ch, cw coarse_t1.shape[:2] change_map np.zeros((h, w, bands)) # 保存每个细像元的变化量 for i in range(ch): for j in range(cw): # 当前粗像元覆盖的细像元范围 r0, r1 i * scale, min((i 1) * scale, h) c0, c1 j * scale, min((j 1) * scale, w) cm_sub class_map[r0:r1, c0:c1] valid cm_sub 0 class_counts np.bincount(cm_sub[valid], minlengthN_CLASS) if class_counts.sum() 0: continue A class_counts / class_counts.sum() # 类别比例 # 粗像元光谱差 delta_y coarse_tk[i, j] - coarse_t1[i, j] # shape: (bands,) # 最小二乘解A delta_E delta_y A np.expand_dims(A, axis0) if A.ndim 1 else A delta_E np.linalg.lstsq(A, delta_y, rcondNone)[0] # shape: (bands,) # 把变化量按类别赋给每个细像元 for b in range(bands): change_map[r0:r1, c0:c1, b] delta_E[b][cm_sub] return change_map实际跑的时候你会发现这层循环非常慢因为Python级别的双层循环在高分影像上完全没有工程效率。我的建议是第一次实现先用这种“直球写法”确保逻辑正确跑通一个小区域验证结果后面想提速再用矢量化或numba重写。瓶颈主要集中在np.linalg.lstsq的反复调用上真实场景可以先求每个coarse像元的类别比例矩阵A的伪逆之后直接做矩阵乘。光谱解混返回的change_map是每个细像元的三维变化量。请注意这一步假设了类别比例在t1到tk之间不变FSDAF原论文也默认在不存在土地覆盖类型突变的前提下使用。如果研究区中间经历过火灾、洪水或者城市化新建需要在预处理阶段把这种突变区域单独掩膜掉否则解混结果会被严重高估或低估。4.4 薄板样条残差分配与最终融合初始预测得到后要把它上采样回粗分辨率与真实MODIS tk比较计算残差。这里的关键是残差是粗网格量必须插值到细网格才能加回初始预测。def thin_plate_residual(coarse_residual, scale): ch, cw coarse_residual.shape # 粗网格中心点坐标用像元索引表示即可 ys, xs np.meshgrid(np.arange(ch), np.arange(cw), indexingij) pts np.stack([xs.ravel() * scale scale // 2, ys.ravel() * scale scale // 2], axis-1) rbf RBFInterpolator( pts, coarse_residual.ravel(), kernelthin_plate_spline, smoothing1e-3, ) h, w ch * scale, cw * scale gx, gy np.meshgrid(np.arange(w), np.arange(h), indexingxy) grid_pts np.stack([gx.ravel(), gy.ravel()], axis-1) return rbf(grid_pts).reshape(h, w)smoothing参数是薄板样条的正则化系数。取0时插值曲面必须严格穿过所有控制点对粗影像上的噪声毫无抵抗力调到1e-3到1e-2之间能在“贴合控制点”和“平滑度”之间取得较好平衡。第一次跑建议用1e-3然后观察残差面的最大值和空间分布如果出现明显的孤立“尖刺”说明平滑度不够增大到1e-2。最终融合的主流程把这些函数串起来fine_t1 read_tif(fine_t1_path)[0].transpose(1, 2, 0) coarse_t1 read_tif(coarse_t1_path)[0].transpose(1, 2, 0) coarse_tk read_tif(coarse_tk_path)[0].transpose(1, 2, 0) class_map, endmember get_class_and_endmember(fine_t1, N_CLASS) change_map unmix_change(...) # 得到每个细像元的变化量 base_pred fine_t1 change_map # 初始预测 # 重采样初始预测到粗分辨率 coarse_pred block_mean(base_pred, scale) residual coarse_tk - coarse_pred residual_fine thin_plate_residual(residual, scale) fusion base_pred residual_fine write_tif(out_path, fusion.transpose(2, 0, 1))block_mean是把细影像按scale窗口取均值得到粗分辨率预测图。要注意MODIS和Landsat的真实像元空间响应函数不一致MODIS的像元并不是严格的矩形box在最高精度实验里应该用MODIS的PSF做卷积降尺度但绝大部分场景下block_mean已经够用。追求更严谨的可以查一下MODIS的LAC任职文档沿用它的band-dependent modulation transfer function。5. 实测效果评估与高频踩坑记录5.1 定量评价指标怎么选没有真实tk时刻Landsat时可以用t2时刻的Landsat做验证把“用t1和tk预测出来的t2影像”和真实t2影像逐像元对比。主流的四个指标指标全称评价重点RMSE均方根误差像素级总体误差越小越好ERGASErreur Relative Globale Adimensionnelle de Synthese综合所有波段和空间分辨率的全局误差SSIM结构相似性纹理和结构的保真度越接近1越好SAM光谱角光谱形状的保真度越小越好我的习惯是RMSE和SSIM必须同时看。RMSE低不代表结构清晰有时过度平滑也会让RMSE下降但SSIM很糟糕。如果SSIM明显偏低大概率是相似像元窗口开太大或者薄板样条的平滑系数过高。5.2 从实际运行里踩出来的五个坑第一个坑是网格未严格对齐。这几乎是所有时空融合代码报错的头号原因。出现“尺寸不匹配”或者融合结果里出现规律的棋盘格条纹时先检查两个输入的transform是否一致。rasterio里可以用src.transform是否完全相等判断。第二个坑是MODIS云残留。我在一片热带研究区第一次跑FSDAF时融合结果里莫名其妙出现了一圈亮斑排查后发现是MODIS影像上残留的小块云被当成了真实反射率变化。后续对所有输入影像都先做QA波段掩膜宁可丢掉那一块区域也不用填补值污染整个解混方程。第三个坑是分类数K的选择。K太小多个异质地物被塞进同一类端元光谱不纯解混结果把不同地物的变化搅在一起K太大噪声被当成独立类别残差分配时产生椒盐噪声。经典经验值是均质农田区4-6类城郊结合部10-15类复杂山地15-20类。可以跑一个简单的K值梯度实验看RMSE变化曲线选拐点。第四个坑是薄板样条插值在大影像上太慢。10000乘10000的网格做RBF插值内存占用和运行时间都会爆炸。我的做法是先对残差影像做分块每块几百个粗像元块与块之间留一定重叠分别插值后再拼接重叠区做线性缝。实测速度能提升一个量级精度损失几乎可以忽略。第五个坑是物候变化跨度过大。FSDAF的默认假设是类别比例不变、端元线性变化但如果是夏到冬这种跨季节预测植被的物候变化经常伴随落叶林叶片掉落、农田收获等“类别属性”改变。这种情况下单纯的光谱解混会系统性低估地表变化幅度建议把预测时段控制在生长季内或者引入时相一致性修正模块。6. 收尾我对FSDAF实现的一句话经验自己在实际跑数据中的体会是FSDAF这类融合算法对数据的挑剔程度远超算法本身。代码逻辑理清楚之后真正决定结果上限的百分之八十在预处理坐标系有没有对齐、云掩膜有没有做干净、两个时相之间是否存在突变。与其花大量时间去调分类数K和薄板样条的平滑系数不如先拿一张小范围测试图把整个pipeline跑通确认每一步输出都在物理合理范围内再放到大图上正式运行。最后再分享一个小技巧第一次实现时不要追求一次性把整景影像全跑完先裁一小块200乘200的Landsat加上对应的一两个MODIS像元做单元测试。这样跑一次只要几秒钟可以快速验证分类、解混、插值三个模块的输出是否合理排查bug的效率会高非常多。等所有模块都符合预期再换成全尺寸数据你踩坑的时间至少能缩短一半。本文还有配套的精品资源点击获取
