如何在R中自动识别锯齿状时间序列上升段的斜率符号？

问题背景 #

存在一个锯齿状时间序列，核心上升段由seq(1,50,1)生成，序列包含较强噪声，无法直接判断斜率方向。全局线性拟合因受下降段干扰，多数情况下无法准确识别上升段的斜率正负。

数据示例 #

set.seed(123) # 设置随机种子保证结果可复现
aa <- c(seq(1,50,1), seq(51,5,-10), seq(1,50,1), seq(51,5,-10),
        seq(1,50,1), seq(51,5,-10), seq(1,50,1), seq(51,5,-10),
        seq(1,50,1), seq(51,5,-10), seq(1,50,1), seq(51,5,-10))
bb <- scale(aa) + 0.5*rnorm(length(aa))

失效的全局拟合方法 #

t <- seq(1, length(aa), 1)
v <- lm(bb ~ t, na.action = na.omit)
v_sign <- v$coefficients[2]
# 输出结果：全局斜率接近0，无法反映上升段的正斜率
cat("全局拟合斜率：", v_sign, "\n")

解决方案 #

1. 滑动窗口局部线性拟合 #

通过固定大小的滑动窗口遍历序列，在每个窗口内拟合线性模型计算斜率，筛选对应上升段的窗口斜率：

library(zoo)

window_size <- 10 # 根据上升段长度（50个点）调整窗口大小
slopes <- rollapply(data = bb, width = window_size, FUN = function(x) {
  t_win <- seq(1, window_size)
  lm(x ~ t_win)$coefficients[2]
}, by = 1, align = "center")

# 可视化斜率分布
plot(slopes, type = "l", xlab = "窗口位置", ylab = "局部斜率", main = "滑动窗口局部斜率")
abline(h = 0, col = "red", lty = 2)

# 统计斜率为正的窗口比例（对应上升段的斜率特征）
positive_slope_ratio <- sum(slopes > 0, na.rm = TRUE) / length(slopes[!is.na(slopes)])
cat("局部正斜率窗口比例：", positive_slope_ratio, "\n")

2. 分段检测后单独拟合上升段 #

先通过差分法结合平滑处理检测上升区间，再对这些区间单独计算斜率：

# 计算差分并做中位数滤波去噪
diff_bb <- diff(bb)
smoothed_diff <- runmed(diff_bb, k = 5)

# 识别连续上升区间
is_rising <- smoothed_diff > 0
rising_segments <- rle(is_rising)
rising_indices <- c()
current_pos <- 1
for (i in seq_along(rising_segments$values)) {
  if (rising_segments$values[i]) {
    start <- current_pos
    end <- current_pos + rising_segments$lengths[i] - 1
    rising_indices <- c(rising_indices, start:end)
  }
  current_pos <- current_pos + rising_segments$lengths[i]
}

# 对上升段单独拟合线性模型
if (length(rising_indices) > 0) {
  t_rising <- t[rising_indices]
  bb_rising <- bb[rising_indices]
  model_rising <- lm(bb_rising ~ t_rising)
  rising_slope_sign <- sign(model_rising$coefficients[2])
  cat("上升段拟合斜率正负：", rising_slope_sign, "\n")
}

3. LOESS局部加权回归估计趋势斜率 #

用LOESS拟合序列趋势，再对趋势曲线求导得到斜率，结合已知的上升段位置统计特征：

# 拟合LOESS趋势，span控制平滑度，需根据数据调整
loess_model <- loess(bb ~ t, span = 0.1)
trend <- predict(loess_model)

# 计算趋势曲线的斜率（差分近似导数）
trend_slope <- diff(trend)

# 基于原始序列结构提取上升段索引
original_rising_indices <- c()
cycle_length <- 50 + 6 # 上升段50点，下降段6点
for (i in seq(1, length(aa), cycle_length)) {
  original_rising_indices <- c(original_rising_indices, i:(i+49))
}
# 提取上升段对应的趋势斜率
rising_trend_slopes <- trend_slope[original_rising_indices[-length(original_rising_indices)]]
# 统计正斜率比例
positive_trend_ratio <- sum(rising_trend_slopes > 0) / length(rising_trend_slopes)
cat("上升段LOESS趋势正斜率比例：", positive_trend_ratio, "\n")

# 可视化原始序列、趋势和斜率
par(mfrow = c(2,1))
plot(t, bb, type = "p", col = "gray", main = "原始序列与LOESS趋势")
lines(t, trend, col = "blue", lwd = 2)
plot(t[-1], trend_slope, type = "l", col = "red", main = "LOESS趋势斜率")
abline(h = 0, lty = 2)

内容的提问来源于stack exchange，提问作者mym123