DeepSeek只是图一乐，真学懂WhittakerFilter还得看我这篇文章。本文将详细介绍差分矩阵，推导WhittakerFilter相关公式，并在GEE中用60行代码实现WhittakerFilter。

Whittaker Filter（惠特克滤波器）是一种用于信号处理和数据分析的平滑技术，主要用于去除噪声并提取信号中的趋势或平滑成分。它由E. T. Whittaker在20世纪初提出，广泛应用于时间序列分析、光谱数据处理、生物信号处理等领域。

核心思想

公式

Whittaker Filter的核心思想是通过最小化一个目标函数来实现信号的平滑。该目标函数由两部分组成：

1. 平滑项：衡量信号平滑度的部分，通常通过二阶差分来量化信号的曲率。

2. 保真项：衡量平滑后信号与原始信号之间差异的部分，确保平滑后的信号不会偏离原始信号太多。

目标函数的形式为：

$$Q=\sum _{i=1}^n{\left(y_i-z_i\right)}^2+\lambda \sum _{i=2}^{n-1}{\left({\Delta }^n{z}_i\right)}^2$$

其中：

• $y_i$是原始信号的第$i$个数据点
• $z_i$是平滑后信号的第$i$个数据点
• $\lambda$是平滑参数，控制平滑程度（$\lambda$越大，信号越平滑）
• ${\Delta }^n{z}_i$是$n$阶差分，用于衡量信号的曲率。

特点

1. 灵活性：通过调整平滑参数 $\lambda$，可以控制平滑程度。$\lambda$ 越大，平滑效果越强；$\lambda$越小，平滑效果越弱。
2. 非参数化：不需要假设信号的分布或模型形式，适用于多种类型的数据。
3. 计算高效：可以通过线性代数方法（如Cholesky分解）高效求解。

应用场景

• 时间序列分析：去除噪声，提取趋势。
• 光谱数据处理：平滑光谱曲线，去除噪声。
• 生物信号处理：如心电图（ECG）或脑电图（EEG）信号的平滑。
• 数据压缩：通过平滑减少数据量。

实现

Whittaker Filter的实现通常涉及以下步骤：

构建差分矩阵 $D$

构建差分矩阵是Whittaker Filter中的一个关键步骤。差分矩阵的作用是通过矩阵运算来计算信号的$n$阶差分，从而量化信号的曲率（即平滑度）。下面详细解释这一过程：

1. 什么是差分？

差分是离散信号处理中常用的操作，用于衡量信号的变化率。具体来说：

• 一阶差分：衡量信号中相邻点之间的变化。
  $$\Delta z_i=z_i-z_{i-1}$$
• 二阶差分：衡量一阶差分的变化，即信号的曲率。
  $${\Delta }^2{z}_i=\Delta z_i-\Delta z_{i-1}=z_i-2z_{i-1}+z_{i-2}$$
  二阶差分可以反映信号的平滑程度。如果二阶差分较小，说明信号在该点附近比较平滑；如果二阶差分较大，说明信号在该点附近变化剧烈。

2. 差分矩阵的作用

差分矩阵$D$是一个线性算子，用于将信号 $z$转换为其二阶差分。具体来说：

• 差分矩阵$D$ 是一个稀疏矩阵，其形式取决于差分的阶数和信号的长度。

• 通过矩阵乘法 $D\cdot z$，可以快速计算信号 $z$ 的二阶差分。

3. 差分矩阵的构建方法

(1) 1阶差分矩阵

对于长度为$m$的信号，1阶差分矩阵$D_1^{\left(m\right)}$是一个$\left(m-1\right)\times m$的稀疏矩阵，形式如下：

$$D_1^{\left(m\right)}=\left[\begin{array}{cccccc}-1& 1& 0& 0& \dots & 0\\ 0& -1& 1& 0& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & 0& -1& 1\end{array}\right]$$

(2) 2阶差分矩阵

2阶差分矩阵$D_2^{\left(m\right)}$是一个$\left(m-2\right)\times m$的稀疏矩阵，可以通过1阶差分矩阵递归计算得到：

$$D_2^{\left(m\right)}=D_1^{\left(m-1\right)}\cdot D_1^{\left(m\right)}=\left[\begin{array}{cccccc}1& -2& 1& 0& \dots & 0\\ 0& 1& -2& 1& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & 1& -2& 1\end{array}\right]$$

其中$D_1^{\left(m-1\right)}$是信号长度为$m-1$的1阶差分矩阵,$D_1^{\left(m\right)}$是信号长度为$m$的1阶差分矩阵。

(3) 3阶差分矩阵

3阶差分矩阵$D_3^{\left(m\right)}$是一个$\left(m-3\right)\times m$的稀疏矩阵，可以通过2阶差分矩阵递归计算得到：

$$D_3^{\left(m\right)}=D_1^{\left(m-2\right)}\cdot D_2^{\left(m\right)}=\left[\begin{array}{ccccccc}1& -3& 3& -1& 0& \dots & 0\\ 0& 1& -3& 3& -1& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & 1& -3& 3& -1\end{array}\right]$$

其中$D_1^{\left(m-2\right)}$是信号长度为$m-2$的1阶差分矩阵,$D_2^{\left(m\right)}$是信号长度为$m$的2阶差分矩阵

(4) n阶差分矩阵

n阶差分矩阵$D_n^{\left(m\right)}$是一个$\left(m-n\right)\times m$的稀疏矩阵，可以通过n-1阶差分矩阵递归计算得到：

$$D_n^{\left(m\right)}=D_1^{\left(m-n+1\right)}\cdot D_{n-1}^{\left(m\right)}$$

其中$D_1^{\left(m-n+1\right)}$是信号长度为$m-n+1$的1阶差分矩阵,$D_{n-1}^{\left(m\right)}$是信号长度为$m$的n-1阶差分矩阵

(5) GEE代码实现

function build_diff_matrix(m, n){
    if (n == 1){
      var rowCount = ee.Number(m);
      var identity_mat = ee.Array.identity(m);
      var left = identity_mat.slice(0,0,rowCount.subtract(1));
      var right = identity_mat.slice(0,1,rowCount);
      var D = left.subtract(right);
    }
    else{
      var D_prev = build_diff_matrix(m, n-1)  //构建(n-1)阶差分矩阵
      var D1 = build_diff_matrix(m-n+1, 1)  //构建一阶差分矩阵
      var D = D1.matrixMultiply(D_prev)
    }
    return D;
};
var m = 5; // 信号长度
var n = 3; // 差分阶数
var D3 = build_diff_matrix(m, n);
print(D3);

推导目标函数

1. 将目标函数转为矩阵形式

为了便于计算，目标函数可以表示为矩阵形式：
$$Q={\Vert Y-Z\Vert }^2+\lambda {\Vert DZ\Vert }^2$$
其中：

• $Y$是原始信号的列向量，
• $Z$是平滑信号的列向量，
• $D$是差分矩阵，用于计算二阶差分。

2. 求目标函数对$Z$的偏导

$$\frac{\partial Q}{\partial Z^{\prime }}=-2\left(Y-Z\right)+2\lambda D^{\prime }DZ$$

3. 令目标函数偏导等于0

因为目标函数为凹函数，因此当偏导数为0时，目标函数到达最小值，此时的 $Z$的即为最能平衡平滑项和保真项的取值。因此目标函数的最小化可以通过求解以下线性方程组来实现：
$$\left(I+\lambda D^{T}D\right)Z=Y$$

求$Z$的最小二乘解

通过求解线性方程组$\left(I+\lambda D^{T}D\right)Z=Y$可以得到平滑信号 $Z$，GEE代码实现如下

var lambda = 10;
var m = imageCollection.size();
var n = 3;
var D = build_diff_matrix(m,n);
var DT = D.transpose();
var I = ee.Array.identity(m);
var I_lambda_DT_D = DT.multiply(lambda).matrixMultiply(D).add(I);
var Y = imageCollection.toArray();
I_lambda_DT_D = ee.Image(I_lambda_DT_D);
var Z = I_lambda_DT_D.matrixSolve(Y);

将$Z$的解转为ee.ImageCollection

求解出的$Z$是类型为array的ee.Image，如何将array图像转为ee.ImageCollection是关键

下面提供一种我自己实现的代码供大家参考

//将array图像转为ee.ImageCollection
var bandNames = imageCollection.first().bandNames();
var newBandNames = bandNames.map(function(band){return ee.String(band).cat("_fit")});
var index_list = ee.List.sequence(0, m.subtract(1));
var image_list = imageCollection.toList(m);
var out_list = index_list.map(function(i){
  var image = image_list.get(i);
  var Z_bands = Z.arraySlice(0, ee.Number(i).byte(), ee.Number(i).add(1).byte()).arrayProject([1]).arrayFlatten([newBandNames]);
  return ee.Image(image).addBands(Z_bands);
});

可视化结果

//可视化
var geometry = /* color: #d63000 */ee.Geometry.Point([89.34334560586538, 32.19430206257235]);
print(ui.Chart.image.series(
  out_col.select(["EVI", "EVI_fit","NDVI", "NDVI_fit"]), geometry, ee.Reducer.mean(), 500)
    .setSeriesNames(["EVI", "EVI_fit","NDVI", "NDVI_fit"])
    .setOptions({
      title: 'whittakerFilter',
      lineWidth: 1,
      pointSize: 3,
}));

参考文献

[1] Eilers P H C. A perfect smoother[J]. Analytical chemistry, 2003, 75(14): 3631-3636.
[2] Gallagher N B. Whittaker Smoother[J]. white paper Eigenvector Research, Inc., www. eigenvector. com, 2018.
[3] Khanal N, Matin M A, Uddin K, et al. A comparison of three temporal smoothing algorithms to improve land cover classification: a case study from NEPAL[J]. Remote Sensing, 2020, 12(18): 2888.

完整代码

完整代码请移步公众号
---------------------------------------------------------------------------点这儿↓↓↓↓----------------------------------------------------------------------------------------------