Hampel滤波器

Hampel滤波器是一种基于中值和中值绝对偏差（MAD）的滤波器，旨在识别和去除时间序列数据中的异常值。相对于传统均值和标准差方法，Hampel滤波器对异常值更具鲁棒性

Hampel滤波器的核心在于中值的计算和MAD的求解。中值代表数据的中间值，而MAD度量了数据点与中值之间的离散程度

中值:

对于包含\(N\)个数据点的数据集\(X\)，中值计算公式如下：

• 如果\(N\)为奇数：中值=\(X_{\frac{N+1}{2}}\)

• 如果\(N\)为偶数：中值=\(\frac{1}{2}(X_{\frac{N}{2}}+X_{\frac{N}{2}+1})\)

中值绝对偏差(MAD)：

MAD是每个数据点\(X_i\)与数据集\(X\)的中值之间绝对差的中值：

MAD=中值(|\(X_i\)-中值(\(X\))|)

其中，\(X_i\)表示每一个数据点

Hampel滤波器异常值判定标准：

如果数据点\(X_i\)被认定为异常值（即其绝对偏差除以MAD超过阈值），那么将其替换为数据集\(X\)的中值；否则，保持数据点的原始值

\[X_i= \begin{cases} 中值(X) \quad \rm{if} \frac{|X_i-中值(X)|}{MAD}>阈值\\ X_i \quad \quad \quad \text{otherwise} \end{cases} \]

# 生成包含异常值的正弦波数据
def generate_data_with_outliers():
    time = np.linspace(0, 10, 100)
    signal = np.sin(time) + np.random.normal(0, 0.1, 100)
    # 添加异常值
    outliers_indices = [20, 40, 60, 80]
    outliers_values = [2.0, -1.9, 2.1, -0.5]
    
    for index, value in zip(outliers_indices, outliers_values):
        signal[index] = value
    return time, signal

def hampel(vals_orig, k=5, v=3):
    """
    使用Hampel滤波器去除时间序列中的异常值。
    参数：
        - vals_orig: numpy数组，原始时间序列数据
        - k: 整数，滤波器窗口的大小（半窗口大小为k/2）
        - v: 浮点数，用于异常值检测的阈值
    返回：
        - vals_filt: numpy数组，经过滤波的时间序列数据
        - outliers_indices: list，异常值的索引列表
    """
    # 创建输入数据的副本
    vals_filt = np.copy(vals_orig)
    outliers_indices = []
    # 定义Hampel滤波器函数
    n = len(vals_orig)
    for i in range(k, n - k):
        # 提取窗口
        window = vals_orig[i - k:i + k + 1]
        # 计算中值和中值绝对偏差（MAD）
        median = np.median(window)
        mad = np.median(np.abs(window - median))
        # 检查当前值是否为异常值
        if np.abs(vals_orig[i] - median) > v * mad:
            # 用中值替换异常值
            vals_filt[i] = median
            # 记录异常值的索引
            outliers_indices.append(i)
    return vals_filt, outliers_indices

    # 生成包含异常值的数据
    time, data = generate_data_with_outliers()    
    # 对数据应用Hampel滤波器
    filtered_data, outliers_indices = hampel(data)
    # 绘制原始数据和滤波后的数据
    plt.figure(figsize=(10, 6))
    plt.plot(time, data, label='origenal data')
    plt.scatter(time[outliers_indices], data[outliers_indices], c='red', 
    	marker='o', label='outliers', s=50)
    plt.plot(time, filtered_data, label='data filtered', linestyle='--', color='red')
    plt.fill_between(time, filtered_data, data, color='red', alpha=0.2, label='outliers region')
    plt.title('outliers removal using hampel identifier')
    plt.grid(True)
    plt.legend()
    plt.show()