基于MATLAB的时空自适应反射率融合模型（STARFM）实现

基于MATLAB的时空自适应反射率融合模型（STARFM）实现

---

一、算法原理与流程

STARFM通过结合高时间分辨率（MODIS）和高空间分辨率（Landsat）数据，生成具有高时空分辨率的合成影像。核心公式为：

其中：

• $L_{t0}$：预测时刻的Landsat反射率
• $M_{tk}$：参考时刻的MODIS反射率
• $W_{ijk}$：空间、时间和光谱权重

实现流程：

1. 数据预处理（配准、重采样）
2. 权重计算（空间/时间/光谱距离）
3. 反射率融合（移动窗口加权平均）
4. 后处理（平滑/去云）

---

二、MATLAB实现代码

1. 数据预处理

%% 读取MODIS和Landsat数据（示例）
modis = readgeoraster('MODIS.tif');  % 500m分辨率
landsat = readgeoraster('Landsat.tif');  % 30m分辨率

%% 数据配准与重采样
modis_resampled = imresize(modis, [30,30], 'bilinear');  % 上采样至30m
landsat_aligned = imwarp(landsat, affine2d([1 0 0; 0 1 0; 0 0 1]));  % 仿射变换对齐

%% 大气校正（以MODIS为例）
modis_surf = modis_resampled ./ (1 + 0.00002*modis_resampled);  % 简化大气校正

2. 权重计算

function W = compute_weights(modis, landsat, t_diff)
    % 空间距离权重（高斯函数）
    [rows, cols] = size(landsat);
    [X, Y] = meshgrid(1:cols, 1:rows);
    spatial_dist = sqrt((X - 0.5).^2 + (Y - 0.5).^2);
    W_space = exp(-(spatial_dist/3).^2);  % 窗口半径3像元
    
    % 时间距离权重
    t_diff_days = t_diff / 86400;  % 转换为天数
    W_time = exp(-t_diff_days^2 / (0.5^2));  % 时间窗口0.5天
    
    % 光谱距离权重（欧氏距离）
    spectral_dist = sqrt(sum((modis - landsat).^2, 3));
    W_spectral = 1 ./ (1 + spectral_dist);
    
    % 综合权重
    W = W_space .* W_time .* W_spectral;
    W = W / sum(W(:));  % 归一化
end

3. 反射率融合

%% 定义预测参数
t_ref = datetime('2020-06-01');  % 参考时刻
t_pred = datetime('2020-06-05');  % 预测时刻
t_diff = t_pred - t_ref;

%% 计算权重矩阵
W = compute_weights(modis_surf, landsat_aligned, t_diff);

%% 移动窗口融合
window_size = 5;  % 5x5窗口
pad = floor(window_size/2);
landsat_pad = padarray(landsat_aligned, [pad pad], 'replicate');

predicted = zeros(size(landsat_aligned));
for i = 1:size(landsat_aligned,1)
    for j = 1:size(landsat_aligned,2)
        % 提取窗口区域
        window = landsat_pad(i:i+2*pad, j:j+2*pad);
        weights = W(i:i+2*pad, j:j+2*pad);
        
        % 加权平均
        predicted(i,j) = sum(sum(window .* weights)) / sum(sum(weights));
    end
end

4. 后处理

%% 空间平滑（非局部均值滤波）
predicted_smooth = nlfilter(predicted, [3 3], @(x) mean(x(:)));

%% 云掩膜处理
qa = readgeoraster('Landsat_QA.tif');
cloud_mask = qa == 0;  % 假设QA值为0表示无云
predicted_final = predicted_smooth .* double(~cloud_mask);

---

三、关键改进

1. 异质区域优化

   引入超分辨率重建（参考文献）：

   % 使用稀疏表示增强MODIS细节
   D = train_dictionary(modis_surf);  % 训练字典
   enhanced_modis = sparse_encode(modis_surf, D);
   

2. 多基图像融合

   结合ESTARFM思想（参考文献）：

   % 使用两对基图像计算权重
   W1 = compute_weights(modis1, landsat1, t_diff);
   W2 = compute_weights(modis2, landsat2, t_diff);
   predicted = (W1*landsat1 + W2*landsat2) / (W1+W2);
   

3. GPU加速

   利用并行计算工具箱：

   gpu_modis = gpuArray(modis_surf);
   gpu_predicted = pagefun(@(x,y,z) compute_weights(x,y,z), gpu_modis, gpu_landsat, t_diff);
   

---

四、验证与可视化

%% 真实影像对比
figure;
subplot(1,3,1); imshow(landsat_aligned); title('真实Landsat');
subplot(1,3,2); imshow(predicted_smooth); title('STARFM预测');
subplot(1,3,3); imshow(imabsdiff(landsat_aligned, predicted_smooth)); title('差异图');

%% 统计指标计算
rmse = sqrt(mean((landsat_aligned(:) - predicted_smooth(:)).^2));
r2 = 1 - sum((landsat_aligned(:)-predicted_smooth(:)).^2) / ...
       sum((landsat_aligned(:)-mean(landsat_aligned(:))).^2);
disp(['RMSE: ', num2str(rmse), ' R²: ', num2str(r2)]);

参考代码 使用Matlab实现遥感时空融合算法STARFM www.youowenfan.com/contentcsq/65759.html

五、应用案例

数据集	输入参数	输出结果
黄河三角洲（2020）	MODIS 500m + Landsat 30m	30m分辨率NDVI时序
亚马逊雨林（2019）	MODIS 1km + GF-1 WFV	16m分辨率地表温度
长江中下游（2021）	MODIS 250m + ZY-3	8m分辨率植被覆盖度

---

六、注意事项

1. 数据配准精度：需保证MODIS和Landsat的空间配准误差<1个像元
2. 时间窗口选择：建议时间间隔不超过5天（参考文献）
3. GPU内存限制：处理大区域时建议分块计算（如1000×1000像元块）
4. 云污染处理：推荐结合QA波段和时序插值（参考文献）

---

七、扩展工具

1. ENVI插件：使用STFMBatchTool批量处理
2. Google Earth Engine：云平台实现大规模计算
3. Python接口：通过rasterio库与MATLAB数据交互