简介：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 米，很多地表细节直接被抹平。这两个数据源经常打架——想要高分辨率就没有高时间频率，想要高时间频率就丢了空间细节。ESTARFM（Enhanced 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 负责读写 GeoTIFF，NumPy 负责数组运算，SciPy 主要用到ndimage和stats模块做距离计算和回归。
我一般用 Python 3.8 到 3.10 跑它，再新的版本容易出现 GDAL 的 wheel 包不好找的问题。如果你用的是 Windows，建议直接用 Conda 装 GDAL，比 pip 省心得多：
conda create -n estarfm python=3.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, geotransform=landsat_t1.geotransform, projection=landsat_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(f"RMSE: {rmse:.4f}") print(f"Correlation: {corr:.4f}") print(f"Bias: {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，取误差最小的那一组。这个调参过程与其说是科学，不如说是玄学——不同地表的反射率方差差异太大，没有放之四海而皆准的值。
说一句建立在血泪上的习惯：我早期用这个算法做一整年时序融合时，跳过了交叉验证直接跑了三十六期数据，结果做完才发现某一景基准期影像有云污染没清理干净，导致那一整个月的融合结果全部作废，返工成本极高。从那以后我每次换数据源或换参数组合，都强制先做一个单日期的交叉验证，再决定是否批量推进。这个习惯听起来简单，但它救过我不止一次。希望帮到你。
本文还有配套的精品资源，点击获取
