本文介紹在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算法處理，我們的資料已經平滑了許多。
至此，大功告成。
歡迎關注：瘋狂學習GIS