近年、リモートセンシングデータを活用した時系列解析は、地球科学や環境モニタリングの分野で注目されています。例えば以前紹介した植生指標であるNDVI（Normalized Difference Vegetation Index）は、その経年変化を追跡することで農業、森林管理、気候変動研究など、幅広い分野に利用されています。
しかし、リモートセンシングデータの解析には独特の課題があります。大気条件の変動、センサーの視野角の変化、雲の存在などの要因により、ノイズが混入し、時系列解析の精度を低下させることがあります。このノイズによって、季節変動や長期トレンドなどの重要な情報を隠されてしまいます。
そこで本記事では、これらの課題に対処するための強力なツールとして、HANTS（Harmonic Analysis of Time Series）を紹介します。HANTS法は、時系列データを傾向成分と周期成分に分解し、ロバストな補間技術を用いてノイズを効果的に除去する手法です。この手法で補正することで、データの周期性を保持しながら、異常値や欠損値を適切に処理できるようになります。
Python 3.11.8, numpy 1.26.3, matplotlib 3.8.2
HANTS法とは？
HANTS（Harmonic ANalysis of Time Series）は、リモートセンシングデータの時系列解析に使用される手法の一つで、時系列データを周期的な成分（調和成分）に分解し、ノイズや異常値を除去しながらデータを再構築する手法です。今回はNDVIでよく利用されるZhou et al., 2015の方法をもとに、この方法の基本的な考え方と動作原理を以下に解説します。
HANTSの基本概念
1. 理論的背景
HANTS法は、フーリエ級数展開の考え方を基礎としています。フーリエ級数展開では、任意の周期関数を正弦波と余弦波の和で表現できるという原理を用います。HANTS法は、この原理を時系列データに適用し、データに含まれる様々な周期成分を抽出します。
2. 基本モデル
HANTS法では、時系列データ $y(t_j)$を以下のようにモデル化します：
$$
\tilde{y}(t_j) = a_0 + \sum_{i=1}^{n_f} \left[ a_i \cos(2 \pi f_i t_j) + b_i \sin(2 \pi f_i t_j) \right]
$$
$$
y(t_j) = \tilde{y}(t_j) + \epsilon(t_j)
$$
ここで、
- ỹ(tj): 再構築された時系列
- y(tj): 元の時系列データ
- a0: データの平均値（ゼロ周波数での係数）
- ai, bi: i番目の調和成分の係数
- fi: i番目の周波数
- tj: 観測時間
- ε(tj): 誤差項（ノイズや異常値を含む）
- nf: 使用する周波数の数
3. HANTS法の特徴
- ロバスト性：異常値に強く、欠損値がある場合でも適用可能
- 周波数の柔軟性：解析したい周波数を自由に選択可能
- ギャップフィリング：欠損データの補間が可能
- 計算効率：大規模なデータセットにも効率的に適用可能
4. パラメータ設定の重要性
HANTS法の性能は、以下のパラメータ設定に依存します：
- 考慮する周波数の数と値
- フィッティングの許容誤差
- 最大反復回数
- 有効データの下限閾値
これらのパラメータを適切に設定することで、特定の応用に最適な結果を得ることができます。
HANTSは、その理論的背景と動作原理により、ノイズの多い時系列データから本質的な周期性を抽出し、高品質な時系列データを再構築することができます。特に、リモートセンシングデータのような不規則なサンプリングや欠損値、異常値を含むデータセットに対して効果的です。
この手法は、植生指数（NDVI）や地表面温度（LST）などの環境データの解析に広く利用されており、季節変動や長期トレンドの抽出に役立ちます。
HANTS法の実装方法
下記のコードをhants.pyという名前でDesktop/labcode/ディレクトリに保存します。
import numpy as np
import matplotlib.pyplot as plt
def hants(time_series, time_points, frequencies, num_frequencies, fit_error_tolerance):
num_samples = len(time_series)
def reconstructed_series(t, intercept, cos_coeffs, sin_coeffs):
result = intercept
for i in range(num_frequencies):
result += cos_coeffs[i] * np.cos(2 * np.pi * frequencies[i] * t) + sin_coeffs[i] * np.sin(2 * np.pi * frequencies[i] * t)
return result
design_matrix = np.zeros((num_samples, 2 * num_frequencies + 1))
design_matrix[:, 0] = 1
for i in range(num_frequencies):
design_matrix[:, 2*i+1] = np.cos(2 * np.pi * frequencies[i] * time_points)
design_matrix[:, 2*i+2] = np.sin(2 * np.pi * frequencies[i] * time_points)
coefficients = np.linalg.lstsq(design_matrix, time_series, rcond=None)[0]
intercept, cos_coeffs, sin_coeffs = coefficients[0], coefficients[1::2], coefficients[2::2]
fitted_series = reconstructed_series(time_points, intercept, cos_coeffs, sin_coeffs)
for _ in range(3): # 3 iterations
residuals = time_series - fitted_series
threshold = np.percentile(np.abs(residuals), 100 * (1 - fit_error_tolerance))
valid_mask = np.abs(residuals) <= threshold
valid_time_series = time_series[valid_mask]
valid_time_points = time_points[valid_mask]
valid_design_matrix = design_matrix[valid_mask]
coefficients = np.linalg.lstsq(valid_design_matrix, valid_time_series, rcond=None)[0]
intercept, cos_coeffs, sin_coeffs = coefficients[0], coefficients[1::2], coefficients[2::2]
fitted_series = reconstructed_series(time_points, intercept, cos_coeffs, sin_coeffs)
return fitted_series, intercept, cos_coeffs, sin_coeffs
# Generate synthetic time series data
np.random.seed(42)
num_days = 365
days = np.arange(num_days)
frequencies = np.array([1/365, 2/365, 3/365, 4/365, 6/365, 12/365])
num_frequencies = len(frequencies)
# Generate the true signal
true_signal = 10
for i, freq in enumerate(frequencies):
amplitude = np.random.uniform(1, 5)
phase = np.random.uniform(0, 2*np.pi)
true_signal += amplitude * np.sin(2 * np.pi * freq * days + phase)
# Add random walk and noise
random_walk = np.cumsum(np.random.normal(0, 0.5, num_days))
true_signal += random_walk
noise = np.random.normal(0, 2, num_days)
noisy_signal = true_signal + noise
# Add outliers
outlier_indices = np.random.choice(num_days, 30, replace=False)
noisy_signal[outlier_indices] += np.random.uniform(10, 20, 30) * np.random.choice([-1, 1], 30)
# Apply HANTS
reconstructed_signal, intercept, cos_coeffs, sin_coeffs = hants(noisy_signal.copy(), days, frequencies, num_frequencies, fit_error_tolerance=0.95)
# Plot results
plt.figure(figsize=(15, 8))
plt.plot(days, true_signal, 'g-', label='True Signal', linewidth=2)
plt.plot(days, noisy_signal, 'b.', alpha=0.5, label='Noisy Data with Outliers')
plt.plot(days, reconstructed_signal, 'r-', label='HANTS Reconstruction', linewidth=2)
plt.legend(fontsize=12)
plt.title('HANTS Algorithm Example with Complex Random Time Series', fontsize=16)
plt.xlabel('Day of Year', fontsize=14)
plt.ylabel('Value', fontsize=14)
plt.grid(True)
plt.tight_layout()
plt.show()
# Print coefficients
print("HANTS Coefficients:")
print(f"Intercept: {intercept:.2f}")
for i in range(num_frequencies):
print(f"Frequency {frequencies[i]:.4f}: cos_coeff = {cos_coeffs[i]:.2f}, sin_coeff = {sin_coeffs[i]:.2f}, "
f"Amplitude = {np.sqrt(cos_coeffs[i]**2 + sin_coeffs[i]**2):.2f}, "
f"Phase = {np.arctan2(sin_coeffs[i], cos_coeffs[i]):.2f}")プログラムを実行する
ターミナルを開き、
$ cd Desktop/labcode/と入力し、ディレクトリを移動します。あとは以下のコマンドを実行するだけです。（ $マークは無視してください）
$ python hants.py実行結果
下記のようなウィンドウが開けば成功です。適当な名前を付けて保存しましょう。
この図の緑の線は、ノイズや外れ値を除いた元の周期信号をプロットしており、青の点はこの信号にランダムノイズと外れ値を付与したもので、実際の観測データを模しています。
HANTSによって再構築された時系列データは赤い線で表されています。このデータを確認すると以下の特徴を示しています：
- 青い点で表される元のノイズの多いデータの中心的な傾向を捉えています。
- 急激な変動や外れ値の影響を受けずに、滑らかな曲線を描いています。
- 緑の線で示される真の信号の全体的な形状をよく再現しています。完全一致ではありませんが、主要な周期的変動を捉えています。
コードの解説
上に書いたソースコードの解説をしていきます。
関数の定義と初期設定
def hants(time_series, time_points, frequencies, num_frequencies, fit_error_tolerance):
num_samples = len(time_series)- この関数は時系列データを入力として受け取り、HANTSアルゴリズムを適用します。
再構築関数の定義
def reconstructed_series(t, intercept, cos_coeffs, sin_coeffs):
result = intercept
for i in range(num_frequencies):
result += cos_coeffs[i] * np.cos(2 * np.pi * frequencies[i] * t) + sin_coeffs[i] * np.sin(2 * np.pi * frequencies[i] * t)
return resultこの内部関数は、フーリエ級数を使って時系列を再構築します。
デザイン行列の作成
design_matrix = np.zeros((num_samples, 2 * num_frequencies + 1))
design_matrix[:, 0] = 1
for i in range(num_frequencies):
design_matrix[:, 2*i+1] = np.cos(2 * np.pi * frequencies[i] * time_points)
design_matrix[:, 2*i+2] = np.sin(2 * np.pi * frequencies[i] * time_points)線形回帰のためのデザイン行列を作成します。各列は異なる周波数のサイン波とコサイン波を表します。
初期フィッティング
coefficients = np.linalg.lstsq(design_matrix, time_series, rcond=None)[0]
intercept, cos_coeffs, sin_coeffs = coefficients[0], coefficients[1::2], coefficients[2::2]
fitted_series = reconstructed_series(time_points, intercept, cos_coeffs, sin_coeffs)最小二乗法を使って初期の係数を求め、それを使って初期のフィッティングを行います。
反復処理
for _ in range(3): # 3 iterations
residuals = time_series - fitted_series
threshold = np.percentile(np.abs(residuals), 100 * (1 - fit_error_tolerance))
valid_mask = np.abs(residuals) <= threshold
valid_time_series = time_series[valid_mask]
valid_time_points = time_points[valid_mask]
valid_design_matrix = design_matrix[valid_mask]
coefficients = np.linalg.lstsq(valid_design_matrix, valid_time_series, rcond=None)[0]
intercept, cos_coeffs, sin_coeffs = coefficients[0], coefficients[1::2], coefficients[2::2]
fitted_series = reconstructed_series(time_points, intercept, cos_coeffs, sin_coeffs)この部分が HANTSのキモです。各反復で以下のステップを実行します
- 残差（元のデータと再構築された系列の差）を計算
- 残差の閾値を設定
- 閾値以下の残差を持つデータポイントのみを選択
- 選択されたデータポイントを使って再度フィッティング
- 新しい係数で系列を再構築
結果の返却
return fitted_series, intercept, cos_coeffs, sin_coeffs最終的な再構築された系列と係数を返します。
ランダムな時系列データの生成
# Generate synthetic time series data
np.random.seed(42)
num_days = 365
days = np.arange(num_days)
frequencies = np.array([1/365, 2/365, 3/365, 4/365, 6/365, 12/365])
num_frequencies = len(frequencies)
# Generate the true signal
true_signal = 10
for i, freq in enumerate(frequencies):
amplitude = np.random.uniform(1, 5)
phase = np.random.uniform(0, 2*np.pi)
true_signal += amplitude * np.sin(2 * np.pi * freq * days + phase)
# Add random walk and noise
random_walk = np.cumsum(np.random.normal(0, 0.5, num_days))
true_signal += random_walk
noise = np.random.normal(0, 2, num_days)
noisy_signal = true_signal + noise
# Add outliers
outlier_indices = np.random.choice(num_days, 30, replace=False)
noisy_signal[outlier_indices] += np.random.uniform(10, 20, 30) * np.random.choice([-1, 1], 30)
- 1年分（365日）のデータポイントを生成します。
- 6つの異なる周波数成分を持つ信号を作成します。
- ランダムウォーク成分を追加して信号をより不規則にします。
- ガウシアンノイズを追加します。
- 30個のランダムな外れ値を追加します。
HANTSアルゴリズムの適用
# Apply HANTS
reconstructed_signal, intercept, cos_coeffs, sin_coeffs = \
hants(noisy_signal.copy(), days, frequencies, num_frequencies, fit_error_tolerance=0.95)- 生成されたノイズの多い時系列データにHANTSアルゴリズムを適用します。
結果のプロット
plt.figure(figsize=(15, 8))
plt.plot(ts, true_signal, 'g-', label='True Signal', linewidth=2)
plt.plot(ts, Y, 'b.', alpha=0.5, label='Noisy Data with Outliers')
plt.plot(ts, Y_hat, 'r-', label='HANTS Reconstruction', linewidth=2)
plt.legend(fontsize=12)
plt.title('HANTS Algorithm Example with Complex Random Time Series', fontsize=16)
plt.xlabel('Day of Year', fontsize=14)
plt.ylabel('Value', fontsize=14)
plt.grid(True)
plt.tight_layout()
plt.show()- 元の時系列データ、ノイズデータ、HANTSによる再構築時系列データを1つのグラフにプロットします。
係数の出力
print("HANTS Coefficients:")
print(f"Intercept: {coef[0]:.2f}")
for i in range(NF):
print(f"Frequency {f[i]:.4f}: Amplitude = {np.sqrt(coef[2*i+1]**2 + coef[2*i+2]**2):.2f}, "
f"Phase = {np.arctan2(coef[2*i+2], coef[2*i+1]):.2f}")- HANTSアルゴリズムによって得られた各周波数成分の振幅と位相を出力します。
最後に
今回はHANTS（Harmonic Analysis of Time Series）アルゴリズムについて、その理論的背景から実装方法、そして実際のデータへの適用例まで詳しく紹介しました。HANTSは一度理解すると、様々な時系列データの解析に簡単かつ効果的に適用できる強力なツールです。
特に衛星画像データに対して、HANTSを適用することで、ノイズや欠損値の影響を軽減し、より信頼性の高い時系列解析が可能になります。ぜひ、この記事で学んだHANTSの手法を、様々な地理的範囲や時間スケールの衛星画像データに適用してみてください。
参考文献
- Zhou, J., Jia, L., & Menenti, M. (2015). Reconstruction of global MODIS NDVI time series: Performance of Harmonic ANalysis of Time Series (HANTS). Remote Sensing of Environment, 163, 217-228.
