简介ESTARFM增强时空自适应反射率融合模型是遥感影像融合中的经典算法尤其适用于复杂异构地表区域的反射率重建。这套以Python语言实现的代码包面向地理信息科学和遥感研究人员旨在通过融合多时相影像解决单一传感器在时间或空间分辨率上的不足可服务于地表覆盖监测、农业估产、环境变化分析等实际任务。压缩包共27个文件大小仅5.18MB结构清晰主要包含Python源码脚本、YAML参数配置文件、TXT日志、PDF算法论文说明、DOCX操作指南以及用于测试的HDR头文件和TM/MODIS多时相遥感数据支持从算法阅读到实际运行的完整流程。其中本次更新版本整合了标准与快速两种融合模式前者保证融合质量后者显著提升大规模数据处理效率测试数据可帮助开发者对比不同参数下的融合效果。这套代码包已有58人学习下载适合中高级遥感开发者用作算法复现、参数调优以及真实数据融合实验的参考与工具。1. ESTARFM 与 Python 实现把逐日 MODIS 和 30 米 Landsat 融合成一套时间序列做长时序地表监测的人都绕不开一个尴尬Landsat 30 米分辨率够用但 16 天重访周期加上云遮挡一年能拿到手的清晰影像往往不到十景MODIS 每天都有分辨率却只有 250 米到 500 米很多地表细节直接被抹平。这两个数据源经常打架——想要高分辨率就没有高时间频率想要高时间频率就丢了空间细节。ESTARFMEnhanced Spatial and Temporal Adaptive Reflectance Fusion Model就是专门解决这个矛盾的它用两个基准期的 Landsat-MODIS 数据对结合预测期的 MODIS 影像在滑动窗口里找相似像元并自适应计算权重最终输出预测期 30 米分辨率反射率。这份 Python 代码实现了完整的 ESTARFM 流程适合做植被物候、农业遥感、地表变化检测的从业者也适合想把这套融合逻辑集成进自己预处理管线的人。2. 算法核心与代码结构从双基准期融合到逐日预测的黑匣子拆解ESTARFM 经常被当成黑匣子——输入几对影像就出结果中间发生了什么没人在意。但用这个代码包之前我建议你先弄明白两件事算法在算什么以及代码是怎么组织的。否则参数一调错出来的图你自己都不敢信。2.1 它在算什么三个输入、一个输出的融合逻辑ESTARFM 的输入固定为三组数据。第一组是基准期 t1 的 Landsat 反射率30 米和同一日期 MODIS 反射率重采样到 30 米网格第二组是基准期 t2 的同类型数据对第三组是预测期 tk 的 MODIS 反射率。输出是 tk 时刻的 30 米模拟 Landsat 反射率。核心思路是「找到变化趋势一致的地表像元加权外推」。因为 Landsat 和 MODIS 在同一时刻观测同一地表两者反射率之间通常存在稳定的线性关系而地表在 t1 到 t2 之间的变化趋势可以用 MODIS 在 tk 相对 t1/t2 的变化量来近似。算法会在每个目标像元周围开一个滑动窗口窗口内逐像元计算与中心像元的光谱距离和空间距离筛选出相似像元然后对每个相似像元求权重——光谱越接近、距离越近权重越高。最后把窗口内相似像元的 MODIS 变化量按权重合成加到中心像元的 Landsat 基准值上。这里有个容易被忽略的细节算法里的「自适应」体现在权重计算上。Gao 等人 2006 年提出原始 STARFM 用的是单一权重函数ESTARFM 在此基础上引入了一个转换系数 V用来描述 MODIS 与 Landsat 反射率之间的线性关系在时间上的稳定性。V 是由两个基准期的 Landsat 和 MODIS 数据对拟合出来的写进代码里就是一组回归系数。也就是说代码不是简单地做差值外推而是先对每个相似像元拟合出 MODIS 到 Landsat 的映射再做时间外推。2.2 代码包结构与主流程先看清文件再动手拿到这份python_estarfm_updated_20211105.zip解开后第一件事不是直接跑而是把目录结构过一遍。我拆过不少遥感算法包这个包的目录组织比较常规核心文件如下表所示文件/目录作用estarfm.py算法主模块包含融合主函数search.py相似像元搜索与筛选weights.py权重计算与归一化regression.py线性回归与转换系数 V 的拟合utils.py影像读写、数组填充、NoData 处理demo/示例数据与运行脚本主流程在estarfm.py的estarfm_fuse()函数里逻辑可以概括为四个动作读取两期基准数据对 → 对每个中心像元搜索周围窗口内的相似像元 → 计算每个相似像元的权重和转换系数 → 合成权重并对预测期 MODIS 做外推。如果你在代码里找「核心中的核心」看weights.py里的权重函数就对了——大部分论文里的公式都收敛在这个文件里。建议你先把demo/下的示例数据跑通一遍再去换自己的数据。示例数据一般已经做好了配准和裁剪跑通只是验证环境没问题真正的坑都在你自己的数据上。2.3 关键函数拆解相似像元、权重和回归分别做了什么逐个说关键函数。search.py里的相似像元搜索直接决定了融合结果的质量。它对窗口内每个像元计算两个距离光谱距离是反射率差值的绝对值经波段求和空间距离就是到中心像元的欧氏距离。筛选时有一个阈值参数只有光谱距离小于阈值的像元才被保留。这个阈值在代码里通常写为threshold默认值一般在 0.05 附近但不同地表类型差异很大——裸地和城区建议调高一些植被区可以保持默认。weights.py里计算的是空间权重和时间权重。空间权重由空间距离决定距离越近权重越大时间权重由 MODIS 在 tk 与 t1/t2 的变化一致性决定。两个权重相乘再归一化得到最终权重。regression.py负责拟合 V 系数。它对每个相似像元取 t1 和 t2 两个时刻的 Landsat-MODIS 反射率对做一元线性回归斜率就是 V。这里代码里有一步容易被忽略回归前要检查两组反射率之间的相关系数相关系数太低说明这个像元在两个基准期之间的变化不规则直接参与融合会给结果引入噪声。部分版本会加一个相关系数阈值过滤如果你拿到的版本没做建议自己补上——这能明显减少融合影像上的斑块噪声。# 伪代码示意 estarfm.py 中的主循环逻辑 for center_y in range(margin, rows - margin): for center_x in range(margin, cols - margin): # 提取窗口内所有像元的反射率 window_landsat_t1 get_window(landsat_t1, center_y, center_x, window_size) window_modis_t1 get_window(modis_t1, center_y, center_x, window_size) window_landsat_t2 get_window(landsat_t2, center_y, center_x, window_size) window_modis_t2 get_window(modis_t2, center_y, center_x, window_size) window_modis_tk get_window(modis_tk, center_y, center_x, window_size) # 筛选相似像元光谱距离小于阈值 similar_idx search_similar_pixels( window_landsat_t1, window_modis_t1, center_y, center_x, threshold ) # 对每个相似像元计算权重空间距离 时间一致性 weights compute_weights(similar_idx, window_size) # 对每个相似像元拟合转换系数 V v_coeffs fit_conversion_coefficient( window_landsat_t1, window_landsat_t2, window_modis_t1, window_modis_t2, similar_idx ) # 加权合成预测期 Landsat 反射率 基准值 加权 MODIS 变化量 * V pred landsat_t1[center_y, center_x] \ np.sum(weights * v_coeffs * (window_modis_tk - window_modis_t1)) result[center_y, center_x] pred这段代码是主循环的骨架。margin是窗口半宽意味着融合结果图像的边缘会有一圈没有输出的区域这是正常现象。window_size决定搜索范围取 30 到 60 之间比较常见窗口太小找不到足够多的相似像元窗口太大又会把不同地类的像元纳进来。调节这两个值对结果影响很大但没有绝对正确的值——Windows 大小和数据分辨率强相关。如果你打开estarfm.py看到的不是这样结构化的函数而是几百行揉在一起的脚本别慌它做的事就是上面四步。动手改代码之前建议先在utils.py里确认影像读出来是数组还是按行流式读取——这直接决定你能处理多大面积的影像。3. 从数据预处理到跑通环境配置与核心参数调优代码包拿到手环境配好了数据集读出来了离看到融合结果还有一段路。这一章的内容顺序和你实际操作的顺序完全一致先准备数据再调参数最后跑并验证。每一步都值得认真对待因为遥感融合这事数据预处理的质量比算法本身更能决定成败。3.1 环境依赖与安装Python 版本、GDAL 和数组运算库先说环境。这份 ESTARFM 的 Python 实现在依赖上并不花哨核心就三个GDAL、NumPy、SciPy。GDAL 负责读写 GeoTIFFNumPy 负责数组运算SciPy 主要用到ndimage和stats模块做距离计算和回归。我一般用 Python 3.8 到 3.10 跑它再新的版本容易出现 GDAL 的 wheel 包不好找的问题。如果你用的是 Windows建议直接用 Conda 装 GDAL比 pip 省心得多conda create -n estarfm python3.9 conda activate estarfm conda install -c conda-forge gdal numpy scipy装完以后做一个快速验证打开一个终端输入python -c from osgeo import gdal; print(gdal.__version__)能打印版本号就说明 GDAL 基础没问题。如果你平时在 VSCode 里写代码记得在.vscode/settings.json里把 Python 解释器指到刚才创建的 Conda 环境不然跑起来会报找不到模块。有一个容易被卡住的点GDAL 读取 GeoTIFF 时返回的波段顺序是 BGR 而不是 RGB如果你的源影像波段顺序不一样融合结果会莫名其妙地偏色。建议在预处理阶段就统一波段顺序代码包里utils.py的read_tiff函数通常已经处理了这一点但处理的是「按文件内顺序读取」不会帮你重排波段。3.2 数据准备流程从原始 MODIS 和 Landsat 到可直接喂进算法的输入数据准备是这个项目里最费时间的一步。ESTARFM 要求三个输入具有完全一致的投影、范围、分辨率和行列数实际上就是把 MODIS 重采样到 Landsat 的网格上去。这里我通常的分步做法是第一步对 Landsat 做大气校正。如果你用的是 Collection 2 Level-2 产品反射率波段已经做好了大气校正可以直接用。但要注意检查 QA 波段把云和云影的像元标记出来后面会用上。第二步对 MODIS 做投影转换和重采样。MODIS 原始数据是正弦投影需要转成 UTM 或 Albers然后用 GDAL 的gdalwarp重采样到 Landsat 影像的四至范围和像元大小。重采样方法用双线性内插就行最邻近会把 MODIS 的块状效应带进来。第三步裁剪到完全一致的范围。真实的操作里最省事的办法是拿一个基准数据通常是 Landsat用它的地理信息直接作为模板把 MODIS 重投影后对齐到这幅模板上然后用gdal_translate -projwin裁剪相同范围最后检查行列数是否一致。这里还要强调云检测的重要性。MODIS 的云掩膜产品MOD35可以直接用把云像元在输入影像里赋成 NoData。如果你不加这一步云污染的 MODIS 像元会带着它的「云反射率特征」参与相似像元筛选和权重计算融合结果里会出现一团一团的白斑——后面避坑章节我会再展开。3.3 核心参数与运行命令窗口大小、阈值和波段配置怎么定环境配好、数据准备完毕接下来就是跑通了。大多数版本的 ESTARFM Python 代码包都会有一个demo.py或者main.py里面定义了输入路径和关键参数。下面给出一个典型的调用方式其中集合了我在实际项目里常用的参数组合from estarfm import estarfm_fuse from utils import read_tiff, write_tiff # 读取输入数据均已预处理为同一投影/范围/分辨率 landsat_t1 read_tiff(data/landsat_t1.tif) # 基准期1 Landsat反射率 modis_t1 read_tiff(data/modis_t1.tif) # 基准期1 MODIS反射率 landsat_t2 read_tiff(data/landsat_t2.tif) # 基准期2 Landsat反射率 modis_t2 read_tiff(data/modis_t2.tif) # 基准期2 MODIS反射率 modis_tk read_tiff(data/modis_tk.tif) # 预测期 MODIS反射率 # 融合参数 params { window_size: 30, # 搜索窗口大小奇数单位像元 threshold: 0.05, # 光谱距离阈值决定相似像元的筛选严格程度 band_indices: [0, 1, 2], # 参与融合的波段索引0-based n_threads: 4 # 多线程加速视机器CPU核心数调整 } # 执行融合 result, quality estarfm_fuse( landsat_t1, modis_t1, landsat_t2, modis_t2, modis_tk, **params ) # 写出结果保留原始地理参考信息 write_tiff(data/fused_landsat_tk.tif, result, geotransformlandsat_t1.geotransform, projectionlandsat_t1.projection)这段代码有几个参数需要你根据自己的数据手动调整。window_size是滑动窗口的边长必须为奇数一般取 21 到 51 之间。如果影像空间分辨率是 30 米窗口取 31 大约对应 930 米的地面范围这个尺度对大部分地表类型是合理的。threshold是相似像元筛选光谱距离阈值取值在 0.01 到 0.1 之间越小越严格。如果你发现输出影像上目标地物边缘很模糊多半是阈值太大把不同地类的像元筛进来了如果大片区域是空的NoData说明阈值太小找不到相似像元。band_indices指定用哪几个波段参与距离计算。如果做 NDVI 时序一般用红波段和近红外两个就够了如果用全色波段融合那就只配一个波段。有一个容易翻车的地方多线程n_threads不是所有版本都支持。如果你用的是旧版代码强行传进去会报 TypeError碰到这种情况直接删掉这个参数单线程跑代价只是慢一些不影响结果。3.4 验证输出合理性的快速检查跑通之后先别急着批量处理花两分钟验证一下结果有没有明显问题。打开输出的 GeoTIFF目视检查三点一是地物边界是否清晰融合影像里道路和农田的边界如果糊成一片说明窗口或阈值设置不合理二是是否有大面积 NoData 区域尤其是影像边缘——这是窗口边界效应导致的正常现象但如果是影像中间出现 NoData多半是相似像元搜索失败三是对照预测期的 MODIS 影像确认输出的空间格局和 MODIS 大致一致否则说明融合权重已经失真。我自己的习惯是用 QGIS 加载输出影像和真实 Landsat 影像做一次 rgb 合成对比这一步能省掉后面大量返工时间。4. 避坑排查融合结果出现条纹和空洞的四个典型原因ESTARFM 这个算法本身不复杂但实际跑起来十个报错里九个出在数据环节还有一个出在参数设置上。这一章把我反复踩过的坑集中列出来每一条都按「现象 → 原因 → 解决」的方式写你遇到类似问题时可以直接对号入座。4.1 融合影像出现整条整条的条纹状噪声现象融合结果的局部区域出现沿轨道方向的平行条纹像被刷子刷过一样目视非常明显但 MODIS 原始影像是干净的。原因这是基准期 Landsat 影像自身有条带噪声Landsat 7 ETM SLC-off 时期的数据尤其多或者 MODIS 重采样后出现了拉丝效应。这类噪声在相似像元筛选中被当作「真实地表反射率」参与了权重计算,融合时被放大。解决在数据预处理阶段就把条带影响的像元标记为 NoData 或者用周边像元插值填补。Landsat 7 的条带可以用landsat7_fix_banding这类脚本先处理再做融合。如果你不想在预处理里引入太多额外步骤最简单的办法是换一景无条带的基准期影像两景基准数据都不要用含条带的。4.2 输出影像中心区域出现大片 NoData现象运行日志没有任何报错但融合结果中间出现成片的 NoData 空洞边缘正常。重新设置threshold后空洞范围改变但不消失。原因这个现象通常是基准期 Landsat 和 MODIS 之间存在系统性云污染导致的。算法在滑动窗口内找不到光谱距离足够小的相似像元于是该像元被标记为 NoData。空洞出现在影像中间而不是边缘说明问题不在窗口边界而在「无相似像元」这一层。解决把输入的四幅影像两期 Landsat 加两期 MODIS分别加载到 GIS 里检查同一位置是否有云或阴影覆盖。如果确认是云污染对 MODIS 做更严格的云掩膜对 Landsat 做云和云影掩膜再重新生成输入。如果影像质量确实没问题试着把threshold从 0.05 调到 0.08但注意这会增加混合像元混入的风险属于两害相权取其轻。4.3 结果影像与预测期 MODIS 空间格局对不上现象融合结果里的地物边界位置和同一时期的 MODIS 影像错位了几个像元看起来像影像没有配准。原因这是输入影像配准不一致的典型症状。MODIS 的几何精度在山区和平坦地区表现不一样与 Landsat 之间可能存在一个像元以上的系统性偏移。算法里相似像元搜索基于空间距离的光谱比较一个像元的偏移就足以让权重算错。解决在预处理阶段做一次配准精化。用 GDAL 的gdal_translate加-a_ullr参数手动纠正偏移或者用gdalwarp -et 0.1做亚像元配准。配准完成后把 MODIS 和 Landsat 叠在一起目视检查一遍——这类问题用眼睛看比用脚本检查更靠谱因为影像偏移在数值上可能只有几十米但误差影响会直接反映在融合结果里。4.4 运行速度极慢甚至内存耗尽崩溃现象影像不到 5000×5000 像元跑了一晚上还没出结果或者直接报 MemoryError 退出。原因ESTARFM 的主循环是逐像元操作的滑动窗口内每个像元都要做相似性搜索和权重计算时间复杂度接近 O(N×W²)其中 N 是影像像元数W 是窗口宽度。如果你把window_size设成 51计算量会成倍增长再加上纯 Python 循环没有优化慢是必然的。解决两个方向。第一缩小窗口比如从 51 改回 31牺牲一些空间连续性换取计算速度。第二如果代码支持把主循环用numba的jit装饰器加速这通常能把运行时间缩短到原来的几十分之一。还有一招是分块处理——把影像切成 1024×1024 的块分别融合块之间保留一定的重叠最后再拼起来。重叠区域取窗口半宽即可这样可以完全抵消边缘效应。5. 结果验证与进阶用法交叉验证、批量处理和多传感器迁移融合结果跑出来了别急着存档。花十五分钟做一个数值验证能让你省下后面一整个月的返工时间。5.1 交叉验证拿真实 Landsat 当天影像做精度评估最有效的验证方法是找一个真实存在 Landsat 影像的日期作为预测期用 ESTARFM 融合出模拟的 Landsat再和真实的 Landsat 逐像元对比。计算均方根误差和相关系数import numpy as np # read_tiff 读取真实 Landsat 和融合模拟结果 real read_tiff(data/landsat_tk_true.tif).astype(np.float64) pred read_tiff(output/fused_landsat_tk.tif).astype(np.float64) # 只统计两幅影像都不是 NoData 的像元 valid (real -9999) (pred -9999) np.isfinite(real) np.isfinite(pred) real_v, pred_v real[valid], pred[valid] # RMSE 和相关系数 rmse np.sqrt(np.mean((real_v - pred_v) ** 2)) corr np.corrcoef(real_v, pred_v)[0, 1] bias np.mean(pred_v - real_v) print(fRMSE: {rmse:.4f}) print(fCorrelation: {corr:.4f}) print(fBias: {bias:.4f})判断标准上很多论文给出的 ESTARFM 融合误差在 0.02 到 0.03 之间反射率绝对值。如果你的 RMSE 明显大于这个范围先排查输入影像配准再排查相似像元筛选参数——通常问题不在算法本身而在数据上。相关系数低于 0.8 就该停下来找原因了这种情况下强行批量处理最后分析出的物候趋势会带上系统性偏差。5.2 批量处理长时序脚本化迭代和中间文件管理验证通过后你就可以放心地批量跑时间序列了。常见做法是写一个外层循环逐期调用融合函数。我习惯把每一期的输出文件命名为fused_YYYYMMDD.tif中间结果统一放在tmp/目录下方便断点续跑和排查。批量处理时别忘了做一件事每次运行前先检查 MODIS 的云覆盖比例。如果某一期 MODIS 云覆盖率超过 30%融合结果基本不可用直接跳过比硬跑更有价值。这个判断可以用 MOD35 云掩膜产品统计云像元占比来实现一行np.mean的事。5.3 参数敏感性和多传感器迁移的边界ESTARFM 并不绑定 Landsat-MODIS 这一对传感器。只要两个传感器满足「高空间频率的粗分辨率数据 低空间频率的细分辨率数据」组合算法逻辑可以直接迁移比如 Sentinel-2 加 MODIS、GF-1 加 MODIS。迁移时唯一必须改的是窗口大小和阈值——因为不同传感器的空间分辨率差异会导致相似像元筛选的物理尺度完全不同。我自己的经验是窗口大小取目标传感器一个像元对应的地面尺度乘以 10 到 20 倍。20 米分辨率的 Sentinel-2 配 500 米 MODIS窗口取 41 左右比较合适。阈值则需要做几次敏感性实验——分别设 0.03、0.05、0.08比较验证期的 RMSE取误差最小的那一组。这个调参过程与其说是科学不如说是玄学——不同地表的反射率方差差异太大没有放之四海而皆准的值。说一句建立在血泪上的习惯我早期用这个算法做一整年时序融合时跳过了交叉验证直接跑了三十六期数据结果做完才发现某一景基准期影像有云污染没清理干净导致那一整个月的融合结果全部作废返工成本极高。从那以后我每次换数据源或换参数组合都强制先做一个单日期的交叉验证再决定是否批量推进。这个习惯听起来简单但它救过我不止一次。希望帮到你。本文还有配套的精品资源点击获取 ESTARFM时空融合模型：Python实现MODIS和Landsat高时空分辨率地表反射率
简介ESTARFM增强时空自适应反射率融合模型是遥感影像融合中的经典算法尤其适用于复杂异构地表区域的反射率重建。这套以Python语言实现的代码包面向地理信息科学和遥感研究人员旨在通过融合多时相影像解决单一传感器在时间或空间分辨率上的不足可服务于地表覆盖监测、农业估产、环境变化分析等实际任务。压缩包共27个文件大小仅5.18MB结构清晰主要包含Python源码脚本、YAML参数配置文件、TXT日志、PDF算法论文说明、DOCX操作指南以及用于测试的HDR头文件和TM/MODIS多时相遥感数据支持从算法阅读到实际运行的完整流程。其中本次更新版本整合了标准与快速两种融合模式前者保证融合质量后者显著提升大规模数据处理效率测试数据可帮助开发者对比不同参数下的融合效果。这套代码包已有58人学习下载适合中高级遥感开发者用作算法复现、参数调优以及真实数据融合实验的参考与工具。1. ESTARFM 与 Python 实现把逐日 MODIS 和 30 米 Landsat 融合成一套时间序列做长时序地表监测的人都绕不开一个尴尬Landsat 30 米分辨率够用但 16 天重访周期加上云遮挡一年能拿到手的清晰影像往往不到十景MODIS 每天都有分辨率却只有 250 米到 500 米很多地表细节直接被抹平。这两个数据源经常打架——想要高分辨率就没有高时间频率想要高时间频率就丢了空间细节。ESTARFMEnhanced Spatial and Temporal Adaptive Reflectance Fusion Model就是专门解决这个矛盾的它用两个基准期的 Landsat-MODIS 数据对结合预测期的 MODIS 影像在滑动窗口里找相似像元并自适应计算权重最终输出预测期 30 米分辨率反射率。这份 Python 代码实现了完整的 ESTARFM 流程适合做植被物候、农业遥感、地表变化检测的从业者也适合想把这套融合逻辑集成进自己预处理管线的人。2. 算法核心与代码结构从双基准期融合到逐日预测的黑匣子拆解ESTARFM 经常被当成黑匣子——输入几对影像就出结果中间发生了什么没人在意。但用这个代码包之前我建议你先弄明白两件事算法在算什么以及代码是怎么组织的。否则参数一调错出来的图你自己都不敢信。2.1 它在算什么三个输入、一个输出的融合逻辑ESTARFM 的输入固定为三组数据。第一组是基准期 t1 的 Landsat 反射率30 米和同一日期 MODIS 反射率重采样到 30 米网格第二组是基准期 t2 的同类型数据对第三组是预测期 tk 的 MODIS 反射率。输出是 tk 时刻的 30 米模拟 Landsat 反射率。核心思路是「找到变化趋势一致的地表像元加权外推」。因为 Landsat 和 MODIS 在同一时刻观测同一地表两者反射率之间通常存在稳定的线性关系而地表在 t1 到 t2 之间的变化趋势可以用 MODIS 在 tk 相对 t1/t2 的变化量来近似。算法会在每个目标像元周围开一个滑动窗口窗口内逐像元计算与中心像元的光谱距离和空间距离筛选出相似像元然后对每个相似像元求权重——光谱越接近、距离越近权重越高。最后把窗口内相似像元的 MODIS 变化量按权重合成加到中心像元的 Landsat 基准值上。这里有个容易被忽略的细节算法里的「自适应」体现在权重计算上。Gao 等人 2006 年提出原始 STARFM 用的是单一权重函数ESTARFM 在此基础上引入了一个转换系数 V用来描述 MODIS 与 Landsat 反射率之间的线性关系在时间上的稳定性。V 是由两个基准期的 Landsat 和 MODIS 数据对拟合出来的写进代码里就是一组回归系数。也就是说代码不是简单地做差值外推而是先对每个相似像元拟合出 MODIS 到 Landsat 的映射再做时间外推。2.2 代码包结构与主流程先看清文件再动手拿到这份python_estarfm_updated_20211105.zip解开后第一件事不是直接跑而是把目录结构过一遍。我拆过不少遥感算法包这个包的目录组织比较常规核心文件如下表所示文件/目录作用estarfm.py算法主模块包含融合主函数search.py相似像元搜索与筛选weights.py权重计算与归一化regression.py线性回归与转换系数 V 的拟合utils.py影像读写、数组填充、NoData 处理demo/示例数据与运行脚本主流程在estarfm.py的estarfm_fuse()函数里逻辑可以概括为四个动作读取两期基准数据对 → 对每个中心像元搜索周围窗口内的相似像元 → 计算每个相似像元的权重和转换系数 → 合成权重并对预测期 MODIS 做外推。如果你在代码里找「核心中的核心」看weights.py里的权重函数就对了——大部分论文里的公式都收敛在这个文件里。建议你先把demo/下的示例数据跑通一遍再去换自己的数据。示例数据一般已经做好了配准和裁剪跑通只是验证环境没问题真正的坑都在你自己的数据上。2.3 关键函数拆解相似像元、权重和回归分别做了什么逐个说关键函数。search.py里的相似像元搜索直接决定了融合结果的质量。它对窗口内每个像元计算两个距离光谱距离是反射率差值的绝对值经波段求和空间距离就是到中心像元的欧氏距离。筛选时有一个阈值参数只有光谱距离小于阈值的像元才被保留。这个阈值在代码里通常写为threshold默认值一般在 0.05 附近但不同地表类型差异很大——裸地和城区建议调高一些植被区可以保持默认。weights.py里计算的是空间权重和时间权重。空间权重由空间距离决定距离越近权重越大时间权重由 MODIS 在 tk 与 t1/t2 的变化一致性决定。两个权重相乘再归一化得到最终权重。regression.py负责拟合 V 系数。它对每个相似像元取 t1 和 t2 两个时刻的 Landsat-MODIS 反射率对做一元线性回归斜率就是 V。这里代码里有一步容易被忽略回归前要检查两组反射率之间的相关系数相关系数太低说明这个像元在两个基准期之间的变化不规则直接参与融合会给结果引入噪声。部分版本会加一个相关系数阈值过滤如果你拿到的版本没做建议自己补上——这能明显减少融合影像上的斑块噪声。# 伪代码示意 estarfm.py 中的主循环逻辑 for center_y in range(margin, rows - margin): for center_x in range(margin, cols - margin): # 提取窗口内所有像元的反射率 window_landsat_t1 get_window(landsat_t1, center_y, center_x, window_size) window_modis_t1 get_window(modis_t1, center_y, center_x, window_size) window_landsat_t2 get_window(landsat_t2, center_y, center_x, window_size) window_modis_t2 get_window(modis_t2, center_y, center_x, window_size) window_modis_tk get_window(modis_tk, center_y, center_x, window_size) # 筛选相似像元光谱距离小于阈值 similar_idx search_similar_pixels( window_landsat_t1, window_modis_t1, center_y, center_x, threshold ) # 对每个相似像元计算权重空间距离 时间一致性 weights compute_weights(similar_idx, window_size) # 对每个相似像元拟合转换系数 V v_coeffs fit_conversion_coefficient( window_landsat_t1, window_landsat_t2, window_modis_t1, window_modis_t2, similar_idx ) # 加权合成预测期 Landsat 反射率 基准值 加权 MODIS 变化量 * V pred landsat_t1[center_y, center_x] \ np.sum(weights * v_coeffs * (window_modis_tk - window_modis_t1)) result[center_y, center_x] pred这段代码是主循环的骨架。margin是窗口半宽意味着融合结果图像的边缘会有一圈没有输出的区域这是正常现象。window_size决定搜索范围取 30 到 60 之间比较常见窗口太小找不到足够多的相似像元窗口太大又会把不同地类的像元纳进来。调节这两个值对结果影响很大但没有绝对正确的值——Windows 大小和数据分辨率强相关。如果你打开estarfm.py看到的不是这样结构化的函数而是几百行揉在一起的脚本别慌它做的事就是上面四步。动手改代码之前建议先在utils.py里确认影像读出来是数组还是按行流式读取——这直接决定你能处理多大面积的影像。3. 从数据预处理到跑通环境配置与核心参数调优代码包拿到手环境配好了数据集读出来了离看到融合结果还有一段路。这一章的内容顺序和你实际操作的顺序完全一致先准备数据再调参数最后跑并验证。每一步都值得认真对待因为遥感融合这事数据预处理的质量比算法本身更能决定成败。3.1 环境依赖与安装Python 版本、GDAL 和数组运算库先说环境。这份 ESTARFM 的 Python 实现在依赖上并不花哨核心就三个GDAL、NumPy、SciPy。GDAL 负责读写 GeoTIFFNumPy 负责数组运算SciPy 主要用到ndimage和stats模块做距离计算和回归。我一般用 Python 3.8 到 3.10 跑它再新的版本容易出现 GDAL 的 wheel 包不好找的问题。如果你用的是 Windows建议直接用 Conda 装 GDAL比 pip 省心得多conda create -n estarfm python3.9 conda activate estarfm conda install -c conda-forge gdal numpy scipy装完以后做一个快速验证打开一个终端输入python -c from osgeo import gdal; print(gdal.__version__)能打印版本号就说明 GDAL 基础没问题。如果你平时在 VSCode 里写代码记得在.vscode/settings.json里把 Python 解释器指到刚才创建的 Conda 环境不然跑起来会报找不到模块。有一个容易被卡住的点GDAL 读取 GeoTIFF 时返回的波段顺序是 BGR 而不是 RGB如果你的源影像波段顺序不一样融合结果会莫名其妙地偏色。建议在预处理阶段就统一波段顺序代码包里utils.py的read_tiff函数通常已经处理了这一点但处理的是「按文件内顺序读取」不会帮你重排波段。3.2 数据准备流程从原始 MODIS 和 Landsat 到可直接喂进算法的输入数据准备是这个项目里最费时间的一步。ESTARFM 要求三个输入具有完全一致的投影、范围、分辨率和行列数实际上就是把 MODIS 重采样到 Landsat 的网格上去。这里我通常的分步做法是第一步对 Landsat 做大气校正。如果你用的是 Collection 2 Level-2 产品反射率波段已经做好了大气校正可以直接用。但要注意检查 QA 波段把云和云影的像元标记出来后面会用上。第二步对 MODIS 做投影转换和重采样。MODIS 原始数据是正弦投影需要转成 UTM 或 Albers然后用 GDAL 的gdalwarp重采样到 Landsat 影像的四至范围和像元大小。重采样方法用双线性内插就行最邻近会把 MODIS 的块状效应带进来。第三步裁剪到完全一致的范围。真实的操作里最省事的办法是拿一个基准数据通常是 Landsat用它的地理信息直接作为模板把 MODIS 重投影后对齐到这幅模板上然后用gdal_translate -projwin裁剪相同范围最后检查行列数是否一致。这里还要强调云检测的重要性。MODIS 的云掩膜产品MOD35可以直接用把云像元在输入影像里赋成 NoData。如果你不加这一步云污染的 MODIS 像元会带着它的「云反射率特征」参与相似像元筛选和权重计算融合结果里会出现一团一团的白斑——后面避坑章节我会再展开。3.3 核心参数与运行命令窗口大小、阈值和波段配置怎么定环境配好、数据准备完毕接下来就是跑通了。大多数版本的 ESTARFM Python 代码包都会有一个demo.py或者main.py里面定义了输入路径和关键参数。下面给出一个典型的调用方式其中集合了我在实际项目里常用的参数组合from estarfm import estarfm_fuse from utils import read_tiff, write_tiff # 读取输入数据均已预处理为同一投影/范围/分辨率 landsat_t1 read_tiff(data/landsat_t1.tif) # 基准期1 Landsat反射率 modis_t1 read_tiff(data/modis_t1.tif) # 基准期1 MODIS反射率 landsat_t2 read_tiff(data/landsat_t2.tif) # 基准期2 Landsat反射率 modis_t2 read_tiff(data/modis_t2.tif) # 基准期2 MODIS反射率 modis_tk read_tiff(data/modis_tk.tif) # 预测期 MODIS反射率 # 融合参数 params { window_size: 30, # 搜索窗口大小奇数单位像元 threshold: 0.05, # 光谱距离阈值决定相似像元的筛选严格程度 band_indices: [0, 1, 2], # 参与融合的波段索引0-based n_threads: 4 # 多线程加速视机器CPU核心数调整 } # 执行融合 result, quality estarfm_fuse( landsat_t1, modis_t1, landsat_t2, modis_t2, modis_tk, **params ) # 写出结果保留原始地理参考信息 write_tiff(data/fused_landsat_tk.tif, result, geotransformlandsat_t1.geotransform, projectionlandsat_t1.projection)这段代码有几个参数需要你根据自己的数据手动调整。window_size是滑动窗口的边长必须为奇数一般取 21 到 51 之间。如果影像空间分辨率是 30 米窗口取 31 大约对应 930 米的地面范围这个尺度对大部分地表类型是合理的。threshold是相似像元筛选光谱距离阈值取值在 0.01 到 0.1 之间越小越严格。如果你发现输出影像上目标地物边缘很模糊多半是阈值太大把不同地类的像元筛进来了如果大片区域是空的NoData说明阈值太小找不到相似像元。band_indices指定用哪几个波段参与距离计算。如果做 NDVI 时序一般用红波段和近红外两个就够了如果用全色波段融合那就只配一个波段。有一个容易翻车的地方多线程n_threads不是所有版本都支持。如果你用的是旧版代码强行传进去会报 TypeError碰到这种情况直接删掉这个参数单线程跑代价只是慢一些不影响结果。3.4 验证输出合理性的快速检查跑通之后先别急着批量处理花两分钟验证一下结果有没有明显问题。打开输出的 GeoTIFF目视检查三点一是地物边界是否清晰融合影像里道路和农田的边界如果糊成一片说明窗口或阈值设置不合理二是是否有大面积 NoData 区域尤其是影像边缘——这是窗口边界效应导致的正常现象但如果是影像中间出现 NoData多半是相似像元搜索失败三是对照预测期的 MODIS 影像确认输出的空间格局和 MODIS 大致一致否则说明融合权重已经失真。我自己的习惯是用 QGIS 加载输出影像和真实 Landsat 影像做一次 rgb 合成对比这一步能省掉后面大量返工时间。4. 避坑排查融合结果出现条纹和空洞的四个典型原因ESTARFM 这个算法本身不复杂但实际跑起来十个报错里九个出在数据环节还有一个出在参数设置上。这一章把我反复踩过的坑集中列出来每一条都按「现象 → 原因 → 解决」的方式写你遇到类似问题时可以直接对号入座。4.1 融合影像出现整条整条的条纹状噪声现象融合结果的局部区域出现沿轨道方向的平行条纹像被刷子刷过一样目视非常明显但 MODIS 原始影像是干净的。原因这是基准期 Landsat 影像自身有条带噪声Landsat 7 ETM SLC-off 时期的数据尤其多或者 MODIS 重采样后出现了拉丝效应。这类噪声在相似像元筛选中被当作「真实地表反射率」参与了权重计算,融合时被放大。解决在数据预处理阶段就把条带影响的像元标记为 NoData 或者用周边像元插值填补。Landsat 7 的条带可以用landsat7_fix_banding这类脚本先处理再做融合。如果你不想在预处理里引入太多额外步骤最简单的办法是换一景无条带的基准期影像两景基准数据都不要用含条带的。4.2 输出影像中心区域出现大片 NoData现象运行日志没有任何报错但融合结果中间出现成片的 NoData 空洞边缘正常。重新设置threshold后空洞范围改变但不消失。原因这个现象通常是基准期 Landsat 和 MODIS 之间存在系统性云污染导致的。算法在滑动窗口内找不到光谱距离足够小的相似像元于是该像元被标记为 NoData。空洞出现在影像中间而不是边缘说明问题不在窗口边界而在「无相似像元」这一层。解决把输入的四幅影像两期 Landsat 加两期 MODIS分别加载到 GIS 里检查同一位置是否有云或阴影覆盖。如果确认是云污染对 MODIS 做更严格的云掩膜对 Landsat 做云和云影掩膜再重新生成输入。如果影像质量确实没问题试着把threshold从 0.05 调到 0.08但注意这会增加混合像元混入的风险属于两害相权取其轻。4.3 结果影像与预测期 MODIS 空间格局对不上现象融合结果里的地物边界位置和同一时期的 MODIS 影像错位了几个像元看起来像影像没有配准。原因这是输入影像配准不一致的典型症状。MODIS 的几何精度在山区和平坦地区表现不一样与 Landsat 之间可能存在一个像元以上的系统性偏移。算法里相似像元搜索基于空间距离的光谱比较一个像元的偏移就足以让权重算错。解决在预处理阶段做一次配准精化。用 GDAL 的gdal_translate加-a_ullr参数手动纠正偏移或者用gdalwarp -et 0.1做亚像元配准。配准完成后把 MODIS 和 Landsat 叠在一起目视检查一遍——这类问题用眼睛看比用脚本检查更靠谱因为影像偏移在数值上可能只有几十米但误差影响会直接反映在融合结果里。4.4 运行速度极慢甚至内存耗尽崩溃现象影像不到 5000×5000 像元跑了一晚上还没出结果或者直接报 MemoryError 退出。原因ESTARFM 的主循环是逐像元操作的滑动窗口内每个像元都要做相似性搜索和权重计算时间复杂度接近 O(N×W²)其中 N 是影像像元数W 是窗口宽度。如果你把window_size设成 51计算量会成倍增长再加上纯 Python 循环没有优化慢是必然的。解决两个方向。第一缩小窗口比如从 51 改回 31牺牲一些空间连续性换取计算速度。第二如果代码支持把主循环用numba的jit装饰器加速这通常能把运行时间缩短到原来的几十分之一。还有一招是分块处理——把影像切成 1024×1024 的块分别融合块之间保留一定的重叠最后再拼起来。重叠区域取窗口半宽即可这样可以完全抵消边缘效应。5. 结果验证与进阶用法交叉验证、批量处理和多传感器迁移融合结果跑出来了别急着存档。花十五分钟做一个数值验证能让你省下后面一整个月的返工时间。5.1 交叉验证拿真实 Landsat 当天影像做精度评估最有效的验证方法是找一个真实存在 Landsat 影像的日期作为预测期用 ESTARFM 融合出模拟的 Landsat再和真实的 Landsat 逐像元对比。计算均方根误差和相关系数import numpy as np # read_tiff 读取真实 Landsat 和融合模拟结果 real read_tiff(data/landsat_tk_true.tif).astype(np.float64) pred read_tiff(output/fused_landsat_tk.tif).astype(np.float64) # 只统计两幅影像都不是 NoData 的像元 valid (real -9999) (pred -9999) np.isfinite(real) np.isfinite(pred) real_v, pred_v real[valid], pred[valid] # RMSE 和相关系数 rmse np.sqrt(np.mean((real_v - pred_v) ** 2)) corr np.corrcoef(real_v, pred_v)[0, 1] bias np.mean(pred_v - real_v) print(fRMSE: {rmse:.4f}) print(fCorrelation: {corr:.4f}) print(fBias: {bias:.4f})判断标准上很多论文给出的 ESTARFM 融合误差在 0.02 到 0.03 之间反射率绝对值。如果你的 RMSE 明显大于这个范围先排查输入影像配准再排查相似像元筛选参数——通常问题不在算法本身而在数据上。相关系数低于 0.8 就该停下来找原因了这种情况下强行批量处理最后分析出的物候趋势会带上系统性偏差。5.2 批量处理长时序脚本化迭代和中间文件管理验证通过后你就可以放心地批量跑时间序列了。常见做法是写一个外层循环逐期调用融合函数。我习惯把每一期的输出文件命名为fused_YYYYMMDD.tif中间结果统一放在tmp/目录下方便断点续跑和排查。批量处理时别忘了做一件事每次运行前先检查 MODIS 的云覆盖比例。如果某一期 MODIS 云覆盖率超过 30%融合结果基本不可用直接跳过比硬跑更有价值。这个判断可以用 MOD35 云掩膜产品统计云像元占比来实现一行np.mean的事。5.3 参数敏感性和多传感器迁移的边界ESTARFM 并不绑定 Landsat-MODIS 这一对传感器。只要两个传感器满足「高空间频率的粗分辨率数据 低空间频率的细分辨率数据」组合算法逻辑可以直接迁移比如 Sentinel-2 加 MODIS、GF-1 加 MODIS。迁移时唯一必须改的是窗口大小和阈值——因为不同传感器的空间分辨率差异会导致相似像元筛选的物理尺度完全不同。我自己的经验是窗口大小取目标传感器一个像元对应的地面尺度乘以 10 到 20 倍。20 米分辨率的 Sentinel-2 配 500 米 MODIS窗口取 41 左右比较合适。阈值则需要做几次敏感性实验——分别设 0.03、0.05、0.08比较验证期的 RMSE取误差最小的那一组。这个调参过程与其说是科学不如说是玄学——不同地表的反射率方差差异太大没有放之四海而皆准的值。说一句建立在血泪上的习惯我早期用这个算法做一整年时序融合时跳过了交叉验证直接跑了三十六期数据结果做完才发现某一景基准期影像有云污染没清理干净导致那一整个月的融合结果全部作废返工成本极高。从那以后我每次换数据源或换参数组合都强制先做一个单日期的交叉验证再决定是否批量推进。这个习惯听起来简单但它救过我不止一次。希望帮到你。本文还有配套的精品资源点击获取