
本文介绍一种鲁棒、免参数的时序阶跃检测方法:通过移动平均去噪、趋势消除和自相关谱分析,自动估计均值突变周期k,并绘制分段恒定均值阶梯图;无需预设模型数或依赖bic等模型选择准则。
本文介绍一种鲁棒、免参数的时序阶跃检测方法:通过移动平均去噪、趋势消除和自相关谱分析,自动估计均值突变周期k,并绘制分段恒定均值阶梯图;无需预设模型数或依赖bic等模型选择准则。
在存在加性高斯噪声且方差未知的条件下,直接对原始时间序列进行BIC驱动的多段均值分割(如分段常数建模)易受局部极小值干扰,且BIC本身需预先设定候选分段数,不适用于周期性阶跃结构的自动发现。相比之下,利用信号内在周期性进行频域分析更为稳健高效——因为本例中均值以固定步长K呈规律性下降(每K步减10),该K值即隐含于信号的自相关函数主峰位置。
以下为完整实现流程:
1. 合成带阶跃均值的时序数据
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq, ifft
from scipy.signal import find_peaks
np.random.seed(2)
n_samples = 180
time = np.arange(n_samples)
# 生成初始均值及周期性下降阶梯
base_mean = np.random.randint(60, 90)
K_true = np.random.randint(10, 40) # 真实周期(待估计)
mean = np.full(n_samples, base_mean)
for i in range(K_true, n_samples, K_true):
mean[i:] = mean[i - K_true] - 10
# 添加噪声(方差未知)
noise = np.random.randn(n_samples) * np.abs(np.random.normal(4, 2))
y = mean + noise
2. 信号预处理:去噪与去趋势
- 移动平均滤波(窗口长度建议为K_true量级,如30–50)抑制高频噪声,保留阶跃轮廓;
- 线性趋势消除避免均值漂移干扰周期检测(此处因均值单调下降,需去除整体斜率):
window = 40 ma = np.convolve(y, np.ones(window), mode='valid') / window # 对齐长度:ma比y短window-1,取中心对齐(可选) ma_full = np.concatenate([np.full(window//2, ma[0]), ma, np.full(window//2, ma[-1])])[:len(y)] # 拟合并减去线性趋势 coeffs = np.polyfit(time, ma_full, 1) trend = coeffs[0] * time + coeffs[1] rm_trend = ma_full - trend
3. 自相关分析定位周期K
计算去趋势后信号的自相关函数(ACF),其首个显著峰值对应的滞后即为近似周期K:
corr = np.correlate(rm_trend, rm_trend, mode='full')
corr = corr[len(corr)//2:] # 取正滞后部分
# 寻找首个显著峰值(排除零滞后)
peaks, _ = find_peaks(corr, height=np.max(corr)*0.3, distance=10)
K_est = peaks[0] if len(peaks) > 0 else int(np.round(len(y)/3))
print(f"Estimated step period K ≈ {K_est} (true: {K_true})")
✅ 关键提示:若ACF峰不明显,可改用功率谱(
|FFT|²)分析——对rm_trend做FFT,取幅度平方,主频率对应f₀ = 1/K,故K ≈ 1/f₀(需注意采样率归一化)。
4. 构建并绘制阶梯均值图
利用估计出的K_est,将时间轴划分为长度为K_est的区间,计算每段均值,生成分段常数阶梯线:
def build_step_mean_curve(y, K):
n_segments = len(y) // K
steps = np.zeros(len(y))
for i in range(n_segments):
start, end = i * K, min((i + 1) * K, len(y))
seg_mean = np.mean(y[start:end])
steps[start:end] = seg_mean
# 填充末段不足K的部分
if n_segments * K <h3>总结与注意事项</h3>
- 优势:该方法不依赖噪声方差先验、无需枚举分段点、抗噪性强,特别适合具有明确周期性阶跃结构的信号;
- 局限性:当K过小(
-
进阶建议:若需更高精度,可在粗估K_est后,在
[K_est−5, K_est+5]邻域内微调,重新计算各K对应的BIC(此时BIC作为验证指标而非主算法),选择BIC最优者。
此流程提供了一条从信号特性出发、物理意义清晰、工程落地简便的阶跃检测路径,远优于盲目套用BIC进行暴力搜索。










