基于HANTS算法的时间序列遥感影像去云重构与动态可视化：实现NDVI值异常值剔除与数据平滑处理

基于HANTS的时间序列重构与动态可视化 HANTS 基于最小二乘法和傅立叶变换通过最小二乘法的迭代拟合去除时序NDVI 值中受云污染影响较大的点借助于傅立叶在时间域和频率域的正反变换实现曲线的分解和重构从而达到时序遥感影像去云重构的目的 适用场景: 异常值剔除、填补缺失值、时间序列平滑、数据压缩 输出重构后平滑的栅格tif MATLAB语言替换文件夹路径即可运行窗外的云层总爱和遥感卫星玩捉迷藏NDVI时间序列上冷不丁冒出的异常值就像汤里的老鼠屎。今儿咱们来盘一盘HANTS这个专治不服的算法手把手用MATLAB调教时序数据。HANTS这货本质上是个数学界的黄金搭档——左边搂着最小二乘法右边牵着傅立叶变换。它的套路是先给异常点打上不乖的标签然后通过反复迭代把那些偏离主旋律的波动给摁下去。就像给NDVI曲线做马杀鸡捏掉云污染造成的肿块最后输出丝滑的栅格tif。基于HANTS的时间序列重构与动态可视化 HANTS 基于最小二乘法和傅立叶变换通过最小二乘法的迭代拟合去除时序NDVI 值中受云污染影响较大的点借助于傅立叶在时间域和频率域的正反变换实现曲线的分解和重构从而达到时序遥感影像去云重构的目的 适用场景: 异常值剔除、填补缺失值、时间序列平滑、数据压缩 输出重构后平滑的栅格tif MATLAB语言替换文件夹路径即可运行先上主菜代码把你们MATLAB的路径换成自己的数据夹input_dir E:/哨兵数据/月度NDVI; output_dir E:/重构结果; tif_list dir(fullfile(input_dir,*.tif)); num_images length(tif_list); dates datetime(2020,1:12,15); % 构造时间序列这段先加载12个月的NDVI数据注意这里用平均日期作为时间轴定位点。接下来才是HANTS的核心魔法for row 1:size(ndvi_stack,1) for col 1:size(ndvi_stack,2) ts squeeze(ndvi_stack(row,col,:)); if sum(~isnan(ts)) 5 % 有效观测少于5次直接跳过 recon_ts(row,col,:) NaN; continue end % 傅里叶基函数配置 harmonics 3; % 3阶谐波 base_freq 2*pi/365; % 年周期 A zeros(num_images, 2*harmonics1); A(:,1) 1; for k1:harmonics A(:,2*k) sin(k*base_freq*datenum(dates)); A(:,2*k1) cos(k*base_freq*datenum(dates)); end % 迭代加权最小二乘 weights ones(num_images,1); for iter 1:10 valid weights 0.1; coeff A(valid,:)\ts(valid); recon A*coeff; residuals abs(ts - recon); outliers residuals 0.2*max(ts); % 残差阈值 weights(outliers) 0; end recon_ts(row,col,:) recon; end end这段双重循环逐个像素做时间序列重构。傅里叶基矩阵A构建了时频转换的桥梁3阶谐波足够捕捉年际变化又不会过拟合。加权迭代时每次把残差大的点权重清零相当于给异常值贴封条。动态可视化才是高潮部分上GIF生成代码fig figure(Position,[100 100 800 600]); for m 1:num_images imagesc(squeeze(recon_ts(:,:,m)),[0 1]); colormap(jet); colorbar; title(datestr(dates(m),yyyy-mm)); drawnow frame getframe(fig); im frame2im(frame); [A,map] rgb2ind(im,256); if m 1 imwrite(A,map,ndvi_animation.gif,gif,LoopCount,0,DelayTime,0.5); else imwrite(A,map,ndvi_animation.gif,gif,WriteMode,append,DelayTime,0.5); end end每帧0.5秒的节奏刚好能让植被指数变化像呼吸般自然。从积雪覆盖到盛夏浓绿12个月的轮回在500毫秒间流转那些被云层遮挡的月份也不再是数据黑洞。跑完这串代码记得检查输出tif的元数据是否完整。遇到沿海区域出现锯齿状异常八成是傅里叶阶数设低了把harmonics参数调到4或5试试。数据量太大卡死把外层循环改成分块处理或者启用MATLAB的parfor并行计算。