简介这是一份基于Python实现的FSDAF遥感影像时空融合算法工程包面向遥感与GIS方向的研究者和开发者用于结合Landsat、MODIS等异源多时相影像提升时空分辨率可服务于土地覆盖变化检测与生态监测等场景。资源共361个文件压缩包约7.72MB以311个py源码为核心覆盖数据读取、预处理、特征提取、融合计算与可视化全流程同时包含exe可执行程序、hdr头文件、yaml/cfg配置、pth权重、txt说明及docx算法指导文档并配备环境激活脚本与示例数据便于直接复现与二次开发。FSDAF的实现涉及重采样、自适应融合策略与精度评估代码按功能模块组织用户可借助Python科学计算生态快速迁移至自定义数据集。目前已有1138人学习下载适合希望深入理解时空融合原理并动手实践的中高级Python开发者。 “遥感影像时空融合”这个方向做长时间序列地表监测的人一定不陌生。Landsat空间分辨率够30米但16天重访周期常常被云遮得七零八落MODIS每天都有数据但500米的空间分辨率对精细应用农田分类、城市热岛、流域水质来说又太粗。FSDAFFlexible Spatiotemporal DAta Fusion就是解决这对“空间精度”与“时间密度”矛盾的经典方法之一由Zhu等人2016年发表在IEEE TGRS上直到现在依然是学术圈做时空融合时最常被cue的基线算法。这篇博文我要聊的是怎么用Python把FSDAF从头实现一遍跑通一条完整的数据链路并分享一些我在写代码、调参过程中踩过的坑。内容适合两类人一是刚接触时空融合、想让算法在本地跑起来的研究生或工程师二是看了原文论文、但不知道代码怎么落地的人。我会尽量把每个模块的“为什么这么写”讲透而不是丢一段黑盒代码让你瞎跑。1. 为什么需要时空融合从Landsat和MODIS的互补关系说起1.1 时空分辨率矛盾的典型场景我拿一个典型的农业应用举例。要监测一块区域冬小麦的关键生育期大概需要10天一条的高质量影像。Landsat单颗星16天重访加上云遮挡一个生长季能拿到三四景无云图像就算运气好。MODIS倒是每天都来但500米混合像元里地块边界、田埂、小沟渠全被抹平种植物候变化根本看不出来。时空融合解决的就是这个矛盾假设在t0时刻有一对Landsat和MODIS影像在t1时刻只有MODIS影像我们想通过模型推断出t1时刻的“Landsat级”高分辨率影像。整个过程不需要额外的数据获取只靠算法把MODIS观测到的时间变化信息传递到Landsat的空间细节上。除了农业物候监测这类方法还常用于洪涝灾害快速评估、森林扰动检测、城市热环境分析。最核心的诉求都一样既要空间细节又要时间密度。没有一种真实卫星能同时满足所以只能靠融合。1.2 FSDAF与STARFM等经典算法的差异时空融合并不是FSDAF开的先河。早在2006年Gao等人就提出STARFM它的思路是在t0时刻的高分辨率影像上搜索相似像元把MODIS观测到的时间变化信号按权重传递过去。这个方法实现简单对均质区域效果不错但对异质区域、破碎地物表现很差——因为相似像元数量有限MODIS的混合信号又不干净容易把变化信息过度平滑掉。FSDAF在STARFM的基础上做了两个关键升级第一个是引入土地覆盖分类先把高分辨率影像分成若干类别确保不同地物在后续建模中各归其位第二个是引入混合像元分解的思想在每个低分辨率像元内部建立线性回归关系把MODIS像元级的变化量“分摊”到更细的高分像元上同时还设计了残差再分配机制来补偿预测误差。这套组合拳下来FSDAF在异质性区域的预测准确性和纹理保持能力普遍优于STARFM。2. FSDAF核心原理拆解三步走的预测逻辑2.1 高分辨率影像分类让不同地物各归其位FSDAF的第一件事是把t0时刻的高分辨率影像做分类。注意这里的分类不是提取具体地物类型而是得到一个土地覆盖类别图它决定了后续哪些像元能一起“讨论变化”。为什么要分类因为一个MODIS像元里很可能是混合地物。假设一个像元内同时存在水体、农田和裸土水体在t0到t1之间基本不变农田正在快速生长裸土可能因为翻耕变暗——如果把三类像元混在一起做回归预测值会把三类变化拉成一个不伦不类的平均谁都没预测准。先分类再建模等于给回归加了一个“同质性约束”只有类别相同的像元才会在逻辑上真正属于同一个变化群体。实际实现时KMeans聚类就够用类别数一般取4到8。KMeans收敛快、对初始点不敏感用随机种子固定即可而且sklearn一行调用不需要自己写ISODATA。当然如果你有现成的土地利用分类产品也可以直接作为输入精度会更好但要注意分类类别不要过多否则后续每类样本量不足。2.2 相似像元搜索与权重计算拿到分类图之后FSDAF要为每个高分辨率像元找“同类邻居”。具体来说在以目标像元为中心的固定窗口比如15×15或25×25内只挑选与中心像元属于同一类别、且光谱差值小于某个阈值经常用一个标准差的像元作为相似像元。这一步的动机很实际后面的回归和残差分配都需要样本但又不能什么样本都要。如果窗口里把不同地物像元都拉进来参与运算回归关系会被异质像元污染如果只靠中心像元自己样本量又不够统计意义太弱。相似像元相当于在“数据量”和“同质性”之间取一个平衡点。每个相似像元会获得两项权重空间距离权重和光谱相似性权重。距离越近权重越大光谱越接近权重越大。这样处理之后离中心像元遥远但光谱碰巧接近的像元不会喧宾夺主符合“空间越近越相关”的地学第一定律。2.3 时间变化预测与残差分配时间变化的预测是FSDAF最核心的部分。原始论文的思路是在每个低分辨率像元内部把低分辨率像元值看成内部高分辨率像元在不同类别上的线性混合通过最小二乘回归求出每个类别在t0到t1之间的变化幅度。有了每个类别的时间变化幅度就能把它应用到对应的所有高分辨率像元上得到初步预测影像。但这个初步预测往往存在系统偏差。原因很多比如分类不够精确、回归模型本身有误差、MODIS与Landsat的观测几何不一致。FSDAF的解决办法是加一步残差分配先把初始预测影像按比例聚合回低分辨率尺度得到模拟的MODIS值与真实t1时刻的MODIS影像相减得到每个低分像元的残差然后把这个残差按相似像元的权重和类别占比分配到窗口内的高分辨率像元上加到初步预测结果上。用一句话概括FSDAF的预测链先分类定结构再回归定幅度最后残差修正偏差。三个环节缺一不可这也是它比STARFM多出来的信息处理深度。3. Python实现FSDAF的整体架构与代码详解3.1 环境准备与数据组织FSDAF的依赖库非常少核心就是numpy、scikit-learn和rasterio。numpy负责矩阵运算scikit-learn提供现成的KMeansrasterio负责读写GeoTIFF。GDAL如果装不上rasterio基本能覆盖常用需求。代码里我不会引入任何深度学习框架纯粹用传统机器学习思路实现方便你后续替换成自己的数据处理逻辑。import numpy as np import rasterio from sklearn.cluster import KMeans数据组织是整个流程里最容易翻车的地方。你需要准备三幅完全配准的影像变量含义典型来源h0t0时刻高分辨率影像形状 (h, w, b)Landsat 8/9 SR产品l0t0时刻低分辨率影像形状 (lh, lw, b)MODIS地表反射率产品l1t1时刻低分辨率影像形状 (lh, lw, b)MODIS地表反射率产品这里有两个关键点。第一波段顺序和波段数量必须对齐。Landsat选哪几个波段MODIS就必须有对应波段的等价物一般常用红光、近红外、短波红外这几个对植被变化敏感的组合。第二空间分辨率比例必须明确。假设Landsat是30米MODIS是500米两者比例大约是16.67比1。实际处理时建议先把两组数据重投影到同一坐标系下重采样成整数倍关系比如让一个MODIS像元正好覆盖16×16个Landsat像元否则后面的块分割索引会错位融合结果会出现明显的条纹噪声。3.2 分类与相似像元模块实现分类模块很直接把高分辨率影像的空间维度展开成像素矩阵然后丢给KMeans即可。def classify_high_res(high_array, n_classes6): high_array: (h, w, b) t0时刻高分辨率影像 n_classes: 土地覆盖类别数 返回: (h, w) 分类标签图 h, w, b high_array.shape data high_array.reshape(-1, b) labels KMeans(n_clustersn_classes, random_state42).fit(data).labels_ return labels.reshape(h, w)这段代码只有四行但它是整个FSDAF的基石。分类质量直接决定后续回归和残差分配的可靠性。实际操作中我会建议把分类结果存成TIFF用GIS软件叠加在原始真彩色影像上看一眼确认水体、农田、城镇这些主要地物没有被明显搞混。相似像元搜索我在这里就不贴完整窗口代码了因为纯Python双层循环太慢真实工程中一般用scipy.ndimage的形态学操作或者numba加速。核心逻辑不难在窗口内筛选类别相同、且光谱距离小于阈值的像元然后同时计算空间距离权重和高光谱距离权重。如果你只是先把流程跑通可以先跳过相似像元搜索退化成“同类像元等权分配”融合效果会略降但代码逻辑清晰很多。3.3 残差分配与最终预测这里我给出一份可直接运行的教学版FSDAF主流程。说明一下这是简化版本省略了窗口搜索和相似像元权重用类别等权代替但完整保留了“分类、初步预测、残差分配”三大环节适合理解算法骨架。def fsdaf_simplified(h0, l0, l1, scale16, n_classes6): 教学版FSDAF h0: (h, w, b) t0高分辨率影像 l0: (lh, lw, b) t0低分辨率影像 l1: (lh, lw, b) t1低分辨率影像 scale: 一个低分像元对应的高分像元边长 h, w, b h0.shape lh, lw, _ l1.shape labels classify_high_res(h0, n_classes) h1 np.zeros_like(h0) # 第一阶段按类别传递时间变化量 for idx in range(lh): for jdx in range(lw): si, sj idx * scale, jdx * scale block0 h0[si:siscale, sj:sjscale] block_labs labels[si:siscale, sj:sjscale] block1 block0.copy() delta l1[idx, jdx] - l0[idx, jdx] # 当前低分像元的变化向量 for k in range(n_classes): mask block_labs k if mask.sum() 0: continue block1[mask] block1[mask] delta h1[si:siscale, sj:sjscale] block1 # 第二阶段残差分配 h1_coarse h1.reshape(lh, scale, lw, scale, b).mean(axis(1, 3)) residual l1 - h1_coarse for idx in range(lh): for jdx in range(lw): si, sj idx * scale, jdx * scale block_labs labels[si:siscale, sj:sjscale] for k in range(n_classes): mask block_labs k if mask.sum() 0: continue frac mask.sum() / (scale * scale) region h1[si:siscale, sj:sjscale] region[mask] region[mask] residual[idx, jdx] * frac return h1写这份代码时有一个numpy陷阱必须提一下。执行region h1[si:siscale, sj:sjscale]得到的region是原数组的视图直接对它做region[mask] value这种复合赋值会让布尔索引产生拷贝导致结果根本没有写回h1。要改成region[mask] region[mask] value才能正确更新。我当初在这里排查了很久才发现希望你不用重复踩坑。整个第一阶段做的事情就是“把MODIS观测到的变化量原样加到对应位置的高分像元上”。教学版里每个类别拿到的变化幅度一样实际上在同一低分像元内不同类别地物的变化幅度大概率不同。这个问题正是FSDAF用回归模型去求解各类别变化幅度的原因。第二阶段则是把所有误差在空间上重新摊还尽量减少预测值与真实MODIS观测之间的系统性偏差。4. 实操中的常见问题与排查经验4.1 窗口大小、分类数等关键参数的调参经验窗口大小和分类数是我在跑FSDAF时调得最频繁的两个参数。窗口太小相似像元数量不足回归方程可能不可解或者严重过拟合窗口太大空间局部性丧失预测结果会变得过度平滑细节纹理全丢。建议从15×15起步观察窗口内各类别的像元数量分布如果某类只有不到十个像元就要扩大窗口。类别数通常取4到8。地物越复杂类别数越要多但也不要盲目增加因为每个低分像元内部的高分像元数是有限的类别太多会导致某些类别在一个块内完全没有样本。一个简单经验先跑一遍KMeans并把聚类结果可视化如果水体被拆成了两种边界接壤的类别说明类别数偏多如果城镇和裸土被归为一类说明类别数偏少。还有一个容易被忽视的问题是光谱阈值。FSDAF在筛选相似像元时用标准差作为阈值如果数据本身噪声大或者大气纠正不彻底标准差会被拉大导致相似像元混入太多异质样本。建议在分类之前先做一次简单的异常值剔除把反射率为负和极高值的像元置为无效值。4.2 边界效应与内存优化的处理边界效应是遥感像元级算法绕不开的问题。窗口滑到影像边缘时会越界常见的做法是反射填充或者直接跳过边缘行。跳过边缘意味着融合结果的外围少一圈像元实用中我会用反射填充然后用掩膜把无效范围外的预测结果裁掉这样既不丢失数据也不会在边缘出现伪纹理。内存方面一幅3000×3000、5波段的Landsat影像以float64存储裸数据大约360MB加上中间变量和分类图很容易超过1GB。处理大区域时我习惯分块计算先把整幅影像分成若干个512×512或1024×1024的块每个块独立调用核心算法最后用子模块的权重重叠区拼接。分块计算还有一个好处方便用multiprocessing做并行多核CPU利用率直接拉满。4.3 精度评估与扩展方向融合结果到底准不准我一般同时看两个层面。一是数值精度计算预测影像与真实Landsat影像之间的相关系数、RMSE和ERGAS指数这是论文里最常见的评价指标二是空间纹理把预测影像和真实影像叠在一起做差分图如果差值在空间上呈随机分布说明误差是随机噪声如果出现明显的块状结构说明分类、配准或者残差分配某个环节有系统性问题。跑通教学版之后想进一步贴近论文效果可以考虑三个扩展方向把等权分配改成窗口相似像元加权分配用最小二乘回归替代“全局delta平均分配”求解每类地物的独立变化幅度加上极值像元选择机制保证局部异质区域的变化不被过度平滑。再往后走就是FSDAF 2.0里引入薄板样条插值以及对深度学习大家常说的DSTF一类的模型了但那些都是以理解FSDAF的这套骨架为前提的。现象常见原因处理方法融合结果出现条纹高低分影像几何未配准重投影后重采样保证整数倍关系预测影像过于平滑分类数过少或窗口过大增加类别数、缩小窗口反射率出现负值残差分配过度对结果clip到物理范围0到1局部区域出现马赛克该区域相似像元不足增大窗口或调整光谱阈值数值误差集中在边界未做边缘填充处理使用反射填充或掩膜裁边最后再分享一个我自己很受用的检查习惯。第一次跑通算法后先别急着用整幅影像调参选一个400×400像素的小块试跑把分类图、初步预测、残差分配图、最终预测分别存成GeoTIFF逐个波段叠在原始影像上对比。很多时空融合的问题不是数学没看懂而是中间某个变量在空间上没有对齐。把中间结果可视化一遍比对着公式猜十个小时都管用。如果你跑出来的结果有“马赛克”或“条纹”十有八九不是算法逻辑写错而是Landsat和MODIS在几何配准、重采样这一步出了问题。本文还有配套的精品资源点击获取
Python实现FSDAF遥感影像时空融合：从原理到代码实战
简介这是一份基于Python实现的FSDAF遥感影像时空融合算法工程包面向遥感与GIS方向的研究者和开发者用于结合Landsat、MODIS等异源多时相影像提升时空分辨率可服务于土地覆盖变化检测与生态监测等场景。资源共361个文件压缩包约7.72MB以311个py源码为核心覆盖数据读取、预处理、特征提取、融合计算与可视化全流程同时包含exe可执行程序、hdr头文件、yaml/cfg配置、pth权重、txt说明及docx算法指导文档并配备环境激活脚本与示例数据便于直接复现与二次开发。FSDAF的实现涉及重采样、自适应融合策略与精度评估代码按功能模块组织用户可借助Python科学计算生态快速迁移至自定义数据集。目前已有1138人学习下载适合希望深入理解时空融合原理并动手实践的中高级Python开发者。 “遥感影像时空融合”这个方向做长时间序列地表监测的人一定不陌生。Landsat空间分辨率够30米但16天重访周期常常被云遮得七零八落MODIS每天都有数据但500米的空间分辨率对精细应用农田分类、城市热岛、流域水质来说又太粗。FSDAFFlexible Spatiotemporal DAta Fusion就是解决这对“空间精度”与“时间密度”矛盾的经典方法之一由Zhu等人2016年发表在IEEE TGRS上直到现在依然是学术圈做时空融合时最常被cue的基线算法。这篇博文我要聊的是怎么用Python把FSDAF从头实现一遍跑通一条完整的数据链路并分享一些我在写代码、调参过程中踩过的坑。内容适合两类人一是刚接触时空融合、想让算法在本地跑起来的研究生或工程师二是看了原文论文、但不知道代码怎么落地的人。我会尽量把每个模块的“为什么这么写”讲透而不是丢一段黑盒代码让你瞎跑。1. 为什么需要时空融合从Landsat和MODIS的互补关系说起1.1 时空分辨率矛盾的典型场景我拿一个典型的农业应用举例。要监测一块区域冬小麦的关键生育期大概需要10天一条的高质量影像。Landsat单颗星16天重访加上云遮挡一个生长季能拿到三四景无云图像就算运气好。MODIS倒是每天都来但500米混合像元里地块边界、田埂、小沟渠全被抹平种植物候变化根本看不出来。时空融合解决的就是这个矛盾假设在t0时刻有一对Landsat和MODIS影像在t1时刻只有MODIS影像我们想通过模型推断出t1时刻的“Landsat级”高分辨率影像。整个过程不需要额外的数据获取只靠算法把MODIS观测到的时间变化信息传递到Landsat的空间细节上。除了农业物候监测这类方法还常用于洪涝灾害快速评估、森林扰动检测、城市热环境分析。最核心的诉求都一样既要空间细节又要时间密度。没有一种真实卫星能同时满足所以只能靠融合。1.2 FSDAF与STARFM等经典算法的差异时空融合并不是FSDAF开的先河。早在2006年Gao等人就提出STARFM它的思路是在t0时刻的高分辨率影像上搜索相似像元把MODIS观测到的时间变化信号按权重传递过去。这个方法实现简单对均质区域效果不错但对异质区域、破碎地物表现很差——因为相似像元数量有限MODIS的混合信号又不干净容易把变化信息过度平滑掉。FSDAF在STARFM的基础上做了两个关键升级第一个是引入土地覆盖分类先把高分辨率影像分成若干类别确保不同地物在后续建模中各归其位第二个是引入混合像元分解的思想在每个低分辨率像元内部建立线性回归关系把MODIS像元级的变化量“分摊”到更细的高分像元上同时还设计了残差再分配机制来补偿预测误差。这套组合拳下来FSDAF在异质性区域的预测准确性和纹理保持能力普遍优于STARFM。2. FSDAF核心原理拆解三步走的预测逻辑2.1 高分辨率影像分类让不同地物各归其位FSDAF的第一件事是把t0时刻的高分辨率影像做分类。注意这里的分类不是提取具体地物类型而是得到一个土地覆盖类别图它决定了后续哪些像元能一起“讨论变化”。为什么要分类因为一个MODIS像元里很可能是混合地物。假设一个像元内同时存在水体、农田和裸土水体在t0到t1之间基本不变农田正在快速生长裸土可能因为翻耕变暗——如果把三类像元混在一起做回归预测值会把三类变化拉成一个不伦不类的平均谁都没预测准。先分类再建模等于给回归加了一个“同质性约束”只有类别相同的像元才会在逻辑上真正属于同一个变化群体。实际实现时KMeans聚类就够用类别数一般取4到8。KMeans收敛快、对初始点不敏感用随机种子固定即可而且sklearn一行调用不需要自己写ISODATA。当然如果你有现成的土地利用分类产品也可以直接作为输入精度会更好但要注意分类类别不要过多否则后续每类样本量不足。2.2 相似像元搜索与权重计算拿到分类图之后FSDAF要为每个高分辨率像元找“同类邻居”。具体来说在以目标像元为中心的固定窗口比如15×15或25×25内只挑选与中心像元属于同一类别、且光谱差值小于某个阈值经常用一个标准差的像元作为相似像元。这一步的动机很实际后面的回归和残差分配都需要样本但又不能什么样本都要。如果窗口里把不同地物像元都拉进来参与运算回归关系会被异质像元污染如果只靠中心像元自己样本量又不够统计意义太弱。相似像元相当于在“数据量”和“同质性”之间取一个平衡点。每个相似像元会获得两项权重空间距离权重和光谱相似性权重。距离越近权重越大光谱越接近权重越大。这样处理之后离中心像元遥远但光谱碰巧接近的像元不会喧宾夺主符合“空间越近越相关”的地学第一定律。2.3 时间变化预测与残差分配时间变化的预测是FSDAF最核心的部分。原始论文的思路是在每个低分辨率像元内部把低分辨率像元值看成内部高分辨率像元在不同类别上的线性混合通过最小二乘回归求出每个类别在t0到t1之间的变化幅度。有了每个类别的时间变化幅度就能把它应用到对应的所有高分辨率像元上得到初步预测影像。但这个初步预测往往存在系统偏差。原因很多比如分类不够精确、回归模型本身有误差、MODIS与Landsat的观测几何不一致。FSDAF的解决办法是加一步残差分配先把初始预测影像按比例聚合回低分辨率尺度得到模拟的MODIS值与真实t1时刻的MODIS影像相减得到每个低分像元的残差然后把这个残差按相似像元的权重和类别占比分配到窗口内的高分辨率像元上加到初步预测结果上。用一句话概括FSDAF的预测链先分类定结构再回归定幅度最后残差修正偏差。三个环节缺一不可这也是它比STARFM多出来的信息处理深度。3. Python实现FSDAF的整体架构与代码详解3.1 环境准备与数据组织FSDAF的依赖库非常少核心就是numpy、scikit-learn和rasterio。numpy负责矩阵运算scikit-learn提供现成的KMeansrasterio负责读写GeoTIFF。GDAL如果装不上rasterio基本能覆盖常用需求。代码里我不会引入任何深度学习框架纯粹用传统机器学习思路实现方便你后续替换成自己的数据处理逻辑。import numpy as np import rasterio from sklearn.cluster import KMeans数据组织是整个流程里最容易翻车的地方。你需要准备三幅完全配准的影像变量含义典型来源h0t0时刻高分辨率影像形状 (h, w, b)Landsat 8/9 SR产品l0t0时刻低分辨率影像形状 (lh, lw, b)MODIS地表反射率产品l1t1时刻低分辨率影像形状 (lh, lw, b)MODIS地表反射率产品这里有两个关键点。第一波段顺序和波段数量必须对齐。Landsat选哪几个波段MODIS就必须有对应波段的等价物一般常用红光、近红外、短波红外这几个对植被变化敏感的组合。第二空间分辨率比例必须明确。假设Landsat是30米MODIS是500米两者比例大约是16.67比1。实际处理时建议先把两组数据重投影到同一坐标系下重采样成整数倍关系比如让一个MODIS像元正好覆盖16×16个Landsat像元否则后面的块分割索引会错位融合结果会出现明显的条纹噪声。3.2 分类与相似像元模块实现分类模块很直接把高分辨率影像的空间维度展开成像素矩阵然后丢给KMeans即可。def classify_high_res(high_array, n_classes6): high_array: (h, w, b) t0时刻高分辨率影像 n_classes: 土地覆盖类别数 返回: (h, w) 分类标签图 h, w, b high_array.shape data high_array.reshape(-1, b) labels KMeans(n_clustersn_classes, random_state42).fit(data).labels_ return labels.reshape(h, w)这段代码只有四行但它是整个FSDAF的基石。分类质量直接决定后续回归和残差分配的可靠性。实际操作中我会建议把分类结果存成TIFF用GIS软件叠加在原始真彩色影像上看一眼确认水体、农田、城镇这些主要地物没有被明显搞混。相似像元搜索我在这里就不贴完整窗口代码了因为纯Python双层循环太慢真实工程中一般用scipy.ndimage的形态学操作或者numba加速。核心逻辑不难在窗口内筛选类别相同、且光谱距离小于阈值的像元然后同时计算空间距离权重和高光谱距离权重。如果你只是先把流程跑通可以先跳过相似像元搜索退化成“同类像元等权分配”融合效果会略降但代码逻辑清晰很多。3.3 残差分配与最终预测这里我给出一份可直接运行的教学版FSDAF主流程。说明一下这是简化版本省略了窗口搜索和相似像元权重用类别等权代替但完整保留了“分类、初步预测、残差分配”三大环节适合理解算法骨架。def fsdaf_simplified(h0, l0, l1, scale16, n_classes6): 教学版FSDAF h0: (h, w, b) t0高分辨率影像 l0: (lh, lw, b) t0低分辨率影像 l1: (lh, lw, b) t1低分辨率影像 scale: 一个低分像元对应的高分像元边长 h, w, b h0.shape lh, lw, _ l1.shape labels classify_high_res(h0, n_classes) h1 np.zeros_like(h0) # 第一阶段按类别传递时间变化量 for idx in range(lh): for jdx in range(lw): si, sj idx * scale, jdx * scale block0 h0[si:siscale, sj:sjscale] block_labs labels[si:siscale, sj:sjscale] block1 block0.copy() delta l1[idx, jdx] - l0[idx, jdx] # 当前低分像元的变化向量 for k in range(n_classes): mask block_labs k if mask.sum() 0: continue block1[mask] block1[mask] delta h1[si:siscale, sj:sjscale] block1 # 第二阶段残差分配 h1_coarse h1.reshape(lh, scale, lw, scale, b).mean(axis(1, 3)) residual l1 - h1_coarse for idx in range(lh): for jdx in range(lw): si, sj idx * scale, jdx * scale block_labs labels[si:siscale, sj:sjscale] for k in range(n_classes): mask block_labs k if mask.sum() 0: continue frac mask.sum() / (scale * scale) region h1[si:siscale, sj:sjscale] region[mask] region[mask] residual[idx, jdx] * frac return h1写这份代码时有一个numpy陷阱必须提一下。执行region h1[si:siscale, sj:sjscale]得到的region是原数组的视图直接对它做region[mask] value这种复合赋值会让布尔索引产生拷贝导致结果根本没有写回h1。要改成region[mask] region[mask] value才能正确更新。我当初在这里排查了很久才发现希望你不用重复踩坑。整个第一阶段做的事情就是“把MODIS观测到的变化量原样加到对应位置的高分像元上”。教学版里每个类别拿到的变化幅度一样实际上在同一低分像元内不同类别地物的变化幅度大概率不同。这个问题正是FSDAF用回归模型去求解各类别变化幅度的原因。第二阶段则是把所有误差在空间上重新摊还尽量减少预测值与真实MODIS观测之间的系统性偏差。4. 实操中的常见问题与排查经验4.1 窗口大小、分类数等关键参数的调参经验窗口大小和分类数是我在跑FSDAF时调得最频繁的两个参数。窗口太小相似像元数量不足回归方程可能不可解或者严重过拟合窗口太大空间局部性丧失预测结果会变得过度平滑细节纹理全丢。建议从15×15起步观察窗口内各类别的像元数量分布如果某类只有不到十个像元就要扩大窗口。类别数通常取4到8。地物越复杂类别数越要多但也不要盲目增加因为每个低分像元内部的高分像元数是有限的类别太多会导致某些类别在一个块内完全没有样本。一个简单经验先跑一遍KMeans并把聚类结果可视化如果水体被拆成了两种边界接壤的类别说明类别数偏多如果城镇和裸土被归为一类说明类别数偏少。还有一个容易被忽视的问题是光谱阈值。FSDAF在筛选相似像元时用标准差作为阈值如果数据本身噪声大或者大气纠正不彻底标准差会被拉大导致相似像元混入太多异质样本。建议在分类之前先做一次简单的异常值剔除把反射率为负和极高值的像元置为无效值。4.2 边界效应与内存优化的处理边界效应是遥感像元级算法绕不开的问题。窗口滑到影像边缘时会越界常见的做法是反射填充或者直接跳过边缘行。跳过边缘意味着融合结果的外围少一圈像元实用中我会用反射填充然后用掩膜把无效范围外的预测结果裁掉这样既不丢失数据也不会在边缘出现伪纹理。内存方面一幅3000×3000、5波段的Landsat影像以float64存储裸数据大约360MB加上中间变量和分类图很容易超过1GB。处理大区域时我习惯分块计算先把整幅影像分成若干个512×512或1024×1024的块每个块独立调用核心算法最后用子模块的权重重叠区拼接。分块计算还有一个好处方便用multiprocessing做并行多核CPU利用率直接拉满。4.3 精度评估与扩展方向融合结果到底准不准我一般同时看两个层面。一是数值精度计算预测影像与真实Landsat影像之间的相关系数、RMSE和ERGAS指数这是论文里最常见的评价指标二是空间纹理把预测影像和真实影像叠在一起做差分图如果差值在空间上呈随机分布说明误差是随机噪声如果出现明显的块状结构说明分类、配准或者残差分配某个环节有系统性问题。跑通教学版之后想进一步贴近论文效果可以考虑三个扩展方向把等权分配改成窗口相似像元加权分配用最小二乘回归替代“全局delta平均分配”求解每类地物的独立变化幅度加上极值像元选择机制保证局部异质区域的变化不被过度平滑。再往后走就是FSDAF 2.0里引入薄板样条插值以及对深度学习大家常说的DSTF一类的模型了但那些都是以理解FSDAF的这套骨架为前提的。现象常见原因处理方法融合结果出现条纹高低分影像几何未配准重投影后重采样保证整数倍关系预测影像过于平滑分类数过少或窗口过大增加类别数、缩小窗口反射率出现负值残差分配过度对结果clip到物理范围0到1局部区域出现马赛克该区域相似像元不足增大窗口或调整光谱阈值数值误差集中在边界未做边缘填充处理使用反射填充或掩膜裁边最后再分享一个我自己很受用的检查习惯。第一次跑通算法后先别急着用整幅影像调参选一个400×400像素的小块试跑把分类图、初步预测、残差分配图、最终预测分别存成GeoTIFF逐个波段叠在原始影像上对比。很多时空融合的问题不是数学没看懂而是中间某个变量在空间上没有对齐。把中间结果可视化一遍比对着公式猜十个小时都管用。如果你跑出来的结果有“马赛克”或“条纹”十有八九不是算法逻辑写错而是Landsat和MODIS在几何配准、重采样这一步出了问题。本文还有配套的精品资源点击获取
