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檔案。具體流程如下：

定義了兩個檔案路徑：
- all_file_path：待處理的.csv檔案所在資料夾路徑；
- output_path：儲存處理後資料的資料夾路徑。
使用dir函數獲取指定資料夾中所有以.csv結尾的檔案。
遍歷每個檔案：
- 構建當前檔案的完整路徑。
- 使用readtable函數讀取.csv檔案資料，並保留列名。
- 獲取需要處理的列索引（2到8列）。
- 遍歷這些列索引：
  - 獲取當前列的名稱和資料。
  - 呼叫HANTS函數對列資料進行處理，得到處理後的資料（儲存在yr中）。
  - 將處理後的資料替換原來的列資料。
- 構建儲存處理後資料的檔名，並使用writetable函數將column_data儲存為.csv檔案。

這裡需要注意，HANTS演演算法的幾個引數，大家就依據自己資料的實際情況來設定即可，具體每一個引數的含義在程式碼包中的HANTS.m檔案內就有介紹。通過如上的程式碼，我們即可實現本文的需求。為了進一步驗證HANTS演演算法是否執行正確，我們可以簡單地繪製一個演演算法處理前後的時間序列資料對比圖，如下圖所示。

可以看到，經過HANTS演演算法處理，我們的資料已經平滑了許多。

至此，大功告成。