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()
