近年、リモートセンシングデータの時系列解析は、地球科学や環境モニタリングの分野で重要性を増しています。代表的な例として、植生指標のNDVI（Normalized Difference Vegetation Index）があります。NDVIの経年変化を追跡することで、農業、森林管理、植生調査など、幅広い分野でのモニターや研究に活用されています。
しかし、リモートセンシングデータの解析には特有の課題が存在します。大気条件の変動、センサーの視野角の変化、雲の影響などにより、データにノイズが混入します。このノイズは季節変動や長期トレンドといった重要な情報を不明瞭にし、時系列解析の精度を著しく低下させる可能性があります。
本記事では、これらの課題を効果的に解決するツールとして、Savitzky-Golayフィルタを紹介します。このフィルタは、時系列データから有意な信号を抽出しながら、ノイズを効率的に除去することができます。
Python 3.11.8, numpy 1.26.3, matplotlib 3.8.2, scipy 1.14.1
Savitzky-Golay フィルタとは？
Savitzky-Golay（SG）フィルタは、データの平滑化とノイズ除去によく用いられるデジタルフィルタの一種です。このフィルタは、最小二乗法による多項式近似に基づいてデータの平滑化を行います。具体的には、データの一定区間（ウィンドウ）内で多項式をフィッティングし、その多項式の中心点での値を用いてノイズを低減します。
SGフィルタは、周波数解析だけでなく、時系列データの平滑化にも広く利用されています。特に、信号の重要な特徴（ピークや谷など）を保持しながらノイズを除去できる点が特徴です。
SGフィルタを利用したNDVIの平滑化
SGフィルタは、植生指数であるNDVIの時系列データにも応用されています。一般的にはChen et al., 2004で提案された次式が元になることが多いです。
$$ Y^r_j = \frac{\sum\limits_{i=-m}^{m} C_i \times Y_{i+j}}{N} $$
$Y^r$は再構成されたNDVI、$C_i$ は平滑化のための畳み込み係数、 $Y$はNDVIの生データです。$N$は平滑化ウィンドウサイズで、$N=2m+1$で表され、$m$は平滑化ウィンドウの半幅です。
パラメータの設定
ここで処理するデータに対して２つのパラメータ（ウィンドウサイズと多項式の次数）を設定する必要があります。一般的に、多項式の次数は2～6の範囲で選択されます。ウィンドウサイズはデータの特性に応じて調整し、最適な結果が得られるようにします。
- ウィンドウサイズの影響:
- 大きすぎる場合: 信号の細かな特徴が失われる可能性があります。
- 小さすぎる場合: ノイズが十分に除去されません。
- 多項式の次数の影響:
- 高すぎる場合: 過剰なフィッティングにより、新たなノイズが生じる可能性があります。
- 低すぎる場合: 信号が過度に平滑化され、重要な特徴が失われる可能性があります。
NDVIデータへのSGフィルタの適用方法
SGフィルタの本質は、スライディングウィンドウ内で多項式を使用して元の信号を部分ごとに最小二乗法でフィッティングすることです。具体的なプロセスは以下の通りです。(Chen et al., 2004: Liu et al., 2016)
多項式の設定:
ウィンドウ内のデータポイントに対して、次数$k$の多項式を仮定します。
$$f_k(i) = \sum_{n=0}^{k} b_ni^n, \quad i∈[-m, m]$$
$b_n$はこの多項式の係数で、 $k$は多項式の次数を表します。
最小二乗法による係数の計算:
多項式の係数$b_n$は以下の最小二乗基準を満たすように計算されます。
$$ \frac{\partial}{\partial b_n}\left[\sum_{i=-m}^{m} (f_k(i)-Y_i)^2\right] = 0 $$
これにより、データに最も適合する多項式が得られます。
畳み込み係数の導出:
多項式のn次微分を$i=0$で評価すると、以下の関係式が得られます。
$$ f_k^n(0) = \sum \limits ^m_{i=-m}C^n_iY_i $$
この式により畳み込み係数が求まります。
平滑化の実施:
$N$で割ることで正規化を行い、平滑化後のNDVI値$Y^r$を計算します。
これらの計算式を参考にNDVIの時系列データにSGフィルタで補正を行ってみましょう。
SGフィルタを利用したNDVIの補正の実装
SGフィルタはPythonのscipyライブラリに用意されているので、今回はそれを利用します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
# 提供されたNDVIデータ（1～12月の半月ごとのデータ）
test_ndvi = np.array([
0.1, 0.15,
0.3, 0.2,
0.3, 0.4,
0.5, 0.5,
0.45, 0.3,
0.6, 0.3,
0.6, 0.7,
0.8, 0.75,
0.7, 0.6,
0.4, 0.55,
0.2, 0.4,
0.2, 0.1
])
# データの準備
ndvi_data = test_ndvi
n = len(ndvi_data)
months = np.linspace(1, 12, n)
# フィルタリングパラメータ
window_size = 5
poly_order = 2
# バリデーション
if window_size > n:
raise ValueError("ウィンドウサイズはデータ長以下である必要があります。")
if window_size % 2 == 0:
raise ValueError("ウィンドウサイズは奇数である必要があります。")
if poly_order >= window_size:
raise ValueError("多項式の次数はウィンドウサイズより小さい必要があります。")
# SGフィルタの適用
ndvi_filtered = savgol_filter(ndvi_data, window_size, poly_order)
# 結果のプロット
plt.figure(figsize=(12, 6))
plt.plot(months, ndvi_data, 'ko-', label='Original NDVI', markersize=6)
plt.plot(months, ndvi_filtered, 'b--', label='SG Filtered NDVI', linewidth=2)
plt.xlabel('Month')
plt.ylabel('NDVI Value')
plt.xticks(ticks=range(1, 13), labels=range(1, 13))
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()実行結果
このコードを実行すると以下のグラフが得られます。今回はお試しのためウィンドウサイズを5、次数を2に設定しましたが、この値は、どれだけの平滑化の精度を求めるかで変更してください。
コードの解説
上に書いたソースコードの解説をしていきます。
ステップ 1: ライブラリのインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filterこの部分では、必要なライブラリをインポートしています。numpyは数値計算用、matplotlibはデータのプロット用、savgol_filterはSavitzky-Golayフィルターの適用に使用されます。
ステップ 2: テストデータの準備
test_ndvi = np.array([
0.1, 0.15,
0.3, 0.2,
0.3, 0.4,
0.5, 0.5,
0.45, 0.3,
0.6, 0.3,
0.6, 0.7,
0.8, 0.75,
0.7, 0.6,
0.4, 0.55,
0.2, 0.4,
0.2, 0.1
])
ndvi_data = test_ndvi
n = len(ndvi_data)
months = np.linspace(1, 12, n)test_ndvi: 半月ごとのNDVI値を格納したNumPy配列。各月につき2つのデータポイントがあります。ndvi_data: 処理対象のNDVIデータ。test_ndviを代入しています。n: データポイントの総数（24）。months: 横軸（x軸）用のデータ。1から12までを等間隔に24分割した配列。np.linspace(1, 12, n)は、1から12までの範囲をn個の等間隔の数値に分割します。- 各データポイントが1年のどの時期に対応するかを示します。
ステップ 3: フィルタリングパラメータの設定
# フィルタリングパラメータ
window_size = 5
poly_order = 2window_size（ウィンドウサイズ）: SGフィルタを適用する際のデータウィンドウのサイズ。- 値は 5 に設定されています。
- ウィンドウサイズは奇数でなければなりません。
poly_order（多項式の次数）: ウィンドウ内のデータにフィットさせる多項式の次数。- 値は 2（2次多項式）に設定されています。
- 多項式の次数はウィンドウサイズより小さい必要があります。
ステップ 4: SGフィルタの適用
# SGフィルタの適用
ndvi_filtered = savgol_filter(ndvi_data, window_size, poly_order)savgol_filter関数を使用して、NDVIデータにSGフィルタを適用します。- 引数:
ndvi_data: フィルタを適用するデータ。window_size: ウィンドウサイズ（5）。poly_order: 多項式の次数（2）。
- 出力:
ndvi_filtered: フィルタリング後のデータが格納されます。
- 引数:
- SGフィルタは、データの局所的な領域に多項式をフィットさせ、その中心点の値を用いてデータを平滑化します。
ステップ 5: 結果のプロット
# 結果のプロット
plt.figure(figsize=(12, 6))
plt.plot(months, ndvi_data, 'ko-', label='Original NDVI', markersize=6)
plt.plot(months, ndvi_filtered, 'b--', label='SG Filtered NDVI', linewidth=2)
plt.xlabel('Month')
plt.ylabel('NDVI Value')
plt.xticks(ticks=range(1, 13), labels=range(1, 13))
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()最後に
Savitzky-Golayフィルタは、ノイズを効果的に除去し、データの品質を向上させる強力な手法です。適切なパラメータ設定とともに使用することで、より正確なデータ解析が可能になります。
この手法は衛星データの解析など、さまざまな分野で広く利用されています。ぜひ活用してみてください。
参考文献
- Chen, Jin, et al. “A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter.” Remote sensing of Environment 91.3-4 (2004): 332-344.
- Liu, Yanping, et al. “Applications of savitzky-golay filter for seismic random noise reduction.” Acta Geophysica 64 (2016): 101-124.
