MATLAB时间序列数据重建与平滑：HANTS滤波

本文介绍在MATLAB中，实现基于HANTS算法（时间序列谐波分析法）的长时间序列数据去噪、重建、填补的详细方法。

  HANTS（Harmonic Analysis of Time Series）是一种用于时间序列分析和插值的算法。它基于谐波分析原理，可以从观测数据中提取出周期性变化的信号成分，并进行数据插值和去噪处理。这一算法的主要思想是将时间序列数据分解为多个不同频率的谐波成分，并通过拟合这些成分来重构原始数据。该算法适用于具有任意周期性的时间序列，可以处理缺失值和异常值，并能够保留原始数据的整体趋势和周期性。

  那么在本文中，我们就介绍一下在MATLAB中，基于我们自己的数据，进行HANTS算法处理的方法。

  首先，由于HANTS算法整体非常精密、复杂，因此我们直接下载一位MATLAB用户撰写好的HANTS算法代码包即可，无需自己手动撰写这一部分的代码。下载方法也很简单，大家进入HANTS算法代码包在MATLAB的官方网站即可。进入网站后，如果大家是第一次使用MATLAB的官方网站，需要注册、登录一下自己的账号；随后，选择屏幕右上角的“Download”选项即可；如下图所示。

  下载完毕后，我们将压缩包解压，即可看到如下图所示的文件列表。

  其中，实现HANTS算法的程序其实就是上图中的前两个文件（也就是ApplyHants.m文件与HANTS.m文件），作者将HANTS算法写成了这两个函数，我们在使用时直接调用这两个函数中的一个即可。其中，第一个函数，也就是ApplyHants.m文件对应的函数，适用于输入数据为多维的情况；而如果我们的数据是一维的，例如常见的对NDVI时序数据、遥感反射率时序数据加以重建，那么就用上图中第二个函数，也就是HANTS.m文件对应的函数即可。

  接下来，我们就可以开始对自己的数据加以HANTS算法处理了。在本文中，我们的需求是这样的：在一个文件夹中，包含有大量的.csv文件，其中每一个文件都具有如下图所示的格式。

  其中，第一行为列名，第一列为时间，后面的列都是不同遥感影像波段反射率的时间序列数据。我们希望，对这一文件夹下所有的.csv文件进行遍历，对其中每一个.csv文件的每一列（除了第一列，因为第一列是表示时间的数据）加以HANTS算法处理。

  明确了具体需求，我们就可以开始撰写代码。前面已经提到了，HANTS算法的代码不用我们自己写，就用下载好函数即可；我们只需要将数据读取、数据预处理、结果保存等部分的代码写好，同时按照自己数据的实际情况，配置一下HANTS算法的各个参数即可。

  本文用到的代码如下。

clear;

ni = 414;
nb = 365 * 8 + 361;
nf = 9;
ts = 1 : 8 : (414 * 8 + 1);
HiLo = "none";
low = 0;
high = 1;
fet = 0.1;
dod = 1;
delta = 0.1;
all_file_path = "E:\01_Reflectivity\99_Model_Training\00_Data\02_Extract_Data\16_8DaysSynthesis_After";
output_path = "E:\01_Reflectivity\99_Model_Training\00_Data\02_Extract_Data\17_HANTS";

files = dir(fullfile(all_file_path, "*.csv"));
for i = 1:numel(files)
    file_path = fullfile(all_file_path, files(i).name);
    column_data = readtable(file_path, "ReadVariableNames",true);
    column_name = column_data.Properties.VariableNames;
    column_index = 2 : 8;
    for j = column_index
        one_column_name = column_name{j};
        one_column_data = column_data.(one_column_name);
        [amp, phi, yr] = HANTS(ni, nb, nf, one_column_data, ts, "none", low, high, fet, dod, delta);
%         [amp_Hi, phi_Hi, yr_Hi] = HANTS(ni, nb, nf, one_column_data, ts, "Hi", low, high, fet, dod, delta);
%         [amp_Lo, phi_Lo, yr_Lo] = HANTS(ni, nb, nf, one_column_data, ts, "Lo", low, high, fet, dod, delta);
        column_data.(one_column_name) = yr;
    end
    save_file_name = fullfile(output_path, files(i).name);
    writetable(column_data, save_file_name, "Delimiter", ",");
end

% plot(one_column_data, "b.-");
% hold on;
% plot(yr, "r.-");
% plot(yr_Hi, "k.-");
% plot(yr_Lo, "g.-");
% legend("Original", "none", "Hi", "Lo");

  其中，这段代码的作用是对每个.csv文件中的指定列数据应用HANTS算法进行处理，并将处理后的数据保存为新的.csv文件。具体流程如下：

1. 定义了两个文件路径：
   • all_file_path：待处理的.csv文件所在文件夹路径；
   • output_path：保存处理后数据的文件夹路径。
2. 使用dir函数获取指定文件夹中所有以.csv结尾的文件。
3. 遍历每个文件：
   • 构建当前文件的完整路径。
   • 使用readtable函数读取.csv文件数据，并保留列名。
   • 获取需要处理的列索引（2到8列）。
   • 遍历这些列索引：
     • 获取当前列的名称和数据。
     • 调用HANTS函数对列数据进行处理，得到处理后的数据（存储在yr中）。
     • 将处理后的数据替换原来的列数据。
   • 构建保存处理后数据的文件名，并使用writetable函数将column_data保存为.csv文件。

  这里需要注意，HANTS算法的几个参数，大家就依据自己数据的实际情况来设置即可，具体每一个参数的含义在代码包中的HANTS.m文件内就有介绍。通过如上的代码，我们即可实现本文的需求。为了进一步验证HANTS算法是否执行正确，我们可以简单地绘制一个算法处理前后的时间序列数据对比图，如下图所示。

  可以看到，经过HANTS算法处理，我们的数据已经平滑了许多。