HANTS学习记录

模型介绍

HANTS模型一般用于处理时间序列的植被指数曲线，获得一个平滑的植被生长曲线，同时消除时间序列中的无效值，较小观测中的云、霾等噪声。

HANTS模型的基础是一组傅里叶级数，所以要从傅里叶级数开始，参考知乎文章（傅里叶系列（一）傅里叶级数的推导）如下式：
$$\begin{gathered} y(t)=\frac{a_0}{2}+a_1\cos (\omega t)+b_1\sin (\omega t) \\ +a_2\cos (2\omega t)+b_2\sin (2\omega t) \\ +\dots \\ =\frac{a_0}{2}+\sum _{i=1}^{\infty }\left[a_i\cos \left(i\omega t\right)+b_i\sin \left(i\omega t\right)\right] \end{gathered}$$
式中的 $a_0$、$a_i$、$b_i$分别为不同频率$f_i$的正余弦波，而$\omega =\frac{2\pi }{T}$也被称为角频率（或者角速度，想象一个圆，沿着圆转一圈需要的时间是T，圆的角度是 $2\pi$，角速度就是 $\frac{2\pi }{T}$）。

我们通常所说的频率$f=\frac{1}{T}$，所以角频率与频率的关系是 $\omega =2\pi f$。
因此，上式也可以写作：
$$y(t)=\frac{a_0}{2}+\sum _{i=1}^{\infty }\left[a_i\cos \left(2\pi ift\right)+b_i\sin \left(2\pi ift\right)\right]$$

我们将上式与论文【Reconstruction of global MODIS NDVI time series: Performance of Harmonic ANalysis of Time Series (HANTS)】中的公式(1)做一个对比：
$$y(t_j)=a_0+\sum _{i=1}^m\left[a_i\cos \left(2\pi f_i{t}_j\right)+b_i\sin \left(2\pi f_i{t}_j\right)\right]$$
虽然有一些字母上的上的差别，但是形式是一样的。上式即HANTS模型。

在拟合植被生长周期时，不会将傅里叶级数无穷展开，不然难以计算。并且展的开级数过多，还会造成过拟合（要知道方波也可以用傅里叶级数表示），导致的问题就是，建立的HANTS模型难以适应具有不同生长周期的地表植被（泛化能力弱），一般展开3~4级即可，这里 $m$ 就表示展开的级数，因此一般设置为3~4。

上式中的$f_i$表示第 $i,i=1,2,\dots ,m$ 级正弦、余弦波的频率，其实$f_i=i\times f$，怎么来理解呢？其实傅里叶级数每一级展开中的正弦、余弦波的频率都不相同，随着级数的增大，频率也越高（表现为正余弦波中波峰波谷的加密），但都是基础频率的整数倍，$f$即基础频率$\frac{1}{T}$。

将模型用矩阵表示，就是用矩阵表示，就是：
$$\mathbf{Y}=\mathbf{X}\mathbf{A}$$
其中：
$$\mathbf{Y}={\left[\begin{array}{ccccc}{y}_1& \dots & y_j& \dots & y_N\end{array}\right]}^T$$
是$N\times 1$的矩阵， $N$表示时间序列中的观测值个数。例如，MODIS有8天、16天合成产品，对应的一年中获取的图像数分别是46、23，即观测值的个数$N=46$ 或 $N=23$。
$$\mathbf{A}={\left[\begin{array}{ccccccccc}{a}_0& a_1& b_1& \dots & a_i& b_i& \dots & a_m& b_m\end{array}\right]}^T$$
为模型中的系数矩阵，是$(2m+1)\times 1$的矩阵
$$\mathbf{X}=\left[\begin{array}{ccccccccc}1& \cos (2\pi f{t}_1)& \sin (2\pi f{t}_1)& \dots & \cos (2\pi if{t}_1)& \sin (2\pi if{t}_1)& \dots & \cos (2\pi mf{t}_1)& \sin (2\pi mf{t}_1)\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& \cos (2\pi f{t}_j)& \sin (2\pi f{t}_j)& \dots & \cos (2\pi if{t}_j)& \sin (2\pi if{t}_j)& \dots & \cos (2\pi mf{t}_j)& \sin (2\pi mf{t}_j)\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮& ⋮\\ 1& \cos (2\pi f{t}_N)& \sin (2\pi f{t}_N)& \dots & \cos (2\pi if{t}_N)& \sin (2\pi if{t}_N)& \dots & \cos (2\pi mf{t}_N)& \sin (2\pi mf{t}_N)\end{array}\right]$$
是$N\times (2m+1)$的矩阵。

矩阵中的$f$在前面说过了，是基础频率，与周期$T$有关。那么，模型中周期$T$怎么定义？一般在应用中，将时间序列中观测值的个数$N$作为$T$。一般来说， $N$是确定的，用户设置了参数 $m$ 后，矩阵 $\mathbf{X}$ 可以通过正弦、余弦函数计算出来。
计算的Python代码如下：

import numpy as np
def get_starter_matrix2(nf, ni):
	'''nf 即 m, ni 即 N'''
    nr = min(2 * nf + 1, ni)
    # 计算谐波中的 cos, sin 值，进而利用线性最小二乘拟合模型系数
    mat = np.zeros((ni, nr))
    mat[:, 0] = 1.0
    ang = 2 * np.pi * np.arange(ni) / ni
    for j in np.arange(1, nf + 1):
        mat[:, 2 * j - 1] = np.cos(j * ang)
        mat[:, 2 * j] = np.sin(j * ang)
    
    return mat

（还有一个问题暂时没想清楚。如果我们在分析的时候，只用一年的数据，以一年作为一个周期，可以像上述那样设置，但是如果用10年、20年的时间序列数据作分析的时候，如何定义这个周期？）

参数求解

理想情况

HANTS 模型的应用，关键就是计算系数$a_0$、$a_i$、$b_i$，即求取系数矩阵$\mathbf{A}$。一般采用线性最小二乘方法，解以下方程组：
$${\mathbf{X}}^T\mathbf{Y}-{\mathbf{X}}^T\mathbf{X}\mathbf{A}=0$$
$\mathbf{X}$ 和 $\mathbf{Y}$ 都是已知的，自然可以求解出$\mathbf{A}$。

实际情况

实际情况是，时间序列曲线中，有一些观测值是无效的，比如被云覆盖时，观测值是不能用的。同时，对于植被而言，植被指数也有一定的范围，例如有的文献中将植被NDVI有效的范围$[low,high]$设置为$[0.2,0.8]$，超出这个范围的值，在计算中也是不能用的。因此在计算系数时，要剔除无效值，即将这些位置的方程剔除。

假设用矩阵$\mathbf{p}$ ($N\times 1$)表示观测值是否有效，矩阵中的元素值只有0和1，0表示无效的。那么第一项应表示为${\mathbf{X}}^T(\mathbf{p}\mathbf{Y})$， 第二项表示为 ${\mathbf{X}}^T\mathbf{p}\mathbf{X}$，以上在进行求解系数A。

需要注意的一个问题是，如果有效观测值的个数小于 $2m+1$，方程没法求解。因此要考虑最大可以剔除的无效值个数，一般剔除的无效值个数不超过 $N-(2m+1)$ 就能满足方程的求解。但是，为了使模型更可靠，要求有效观测值的个数要比$2m+1$多出$dod$个，那么，可以剔除的无效值个数就是$N-(2m+1)-dod$。文献【Reconstructing cloudfree NDVI composites using Fourier analysis of time series（这也是提出HANTS的论文）】对$dod$的定义是 Degree of overdeterminedness (DOD): In order to get a more reliable fit the user can decide here to use more data points than the necessary minimum. The minimum number of extra data points, which have to be used in the ultimate fit, is given by the DOD value.