如何在气压数据向量中自动识别显著局部趋势与多段划分

问题背景 #

• 拥有毫巴单位的气压数据向量，常规做法是将过去3小时气压趋势划分为两段，对应9种分类，但实际场景中常存在3段及以上的变化模式（如上升→平稳→大幅下降）
• 已通过find_peaks函数生成包含局部最大值和最小值的maxmin_df，需基于此实现符合阈值（±0.1毫巴/小时）的多段趋势划分

核心思路 #

基于已识别的局部极值点，结合趋势阈值验证，筛选出具有统计显著性的分段点：

1. 补充序列首尾端点到极值列表，确保覆盖完整数据范围
2. 遍历相邻极值点区间，计算区间内的线性趋势斜率并转换为毫巴/小时的变化率
3. 根据阈值判断区间趋势类型（上升/平稳/下降），合并连续同趋势的区间
4. 输出最终的分段结果及对应趋势

实现代码 #

data <- c(966.6003, 966.5772, 966.6102, 966.5933, 966.5768, 966.565, 966.5139, 
966.4826, 966.4626, 966.4595, 966.5143, 966.5732, 966.636, 966.6833, 966.7229, 
966.7517, 966.7523, 966.7955, 966.8255, 966.7964, 966.8406, 966.9269, 966.9753, 
966.9807, 967.001, 967.0519, 967.0389, 967.107, 967.098, 967.0935, 967.0296, 967.0562, 
966.9952, 967.0031, 966.9976, 966.9888, 966.9832, 967.0474, 966.991, 967.0148, 
967.0944, 967.0556, 967.0732, 967.051, 967.1252, 967.0908, 967.0746, 967.0633, 
967.0821, 967.1274, 967.1229, 967.0838, 967.045, 967.0626, 967.0178, 967.024, 967.0806, 
967.0135, 967.0191, 967.0537, 967.0317, 967.0933, 967.0366, 967.0532, 967.0474, 
967.0831, 967.0982, 967.1099, 967.0803, 967.0351, 967.0497, 967.0833, 967.0707, 
967.1096, 967.1709, 967.1645, 967.1979, 967.1406, 967.1523, 967.1288, 967.1345, 
967.1495, 967.1555, 967.1013, 967.0684, 967.1175, 967.1157, 967.0471, 967.0621, 
967.0117, 966.9926, 966.9926, 966.9672, 966.9619, 966.982, 966.9976, 966.9993, 
967.001, 966.9885, 966.9634, 966.9919, 967.0289, 967.0289, 967.037, 967.0547, 967.0158, 
967.0444, 966.9938, 966.9659, 966.9753, 966.9716, 966.9322, 966.9338, 966.8095, 
966.7881, 966.7939, 966.7171, 966.738, 966.7165, 966.7402, 966.7237, 966.6966, 
966.7085, 966.7268, 966.7099, 966.687, 966.6978, 966.6718, 966.6927, 966.7132, 
966.6773, 966.5689, 966.5837, 966.4875, 966.4287, 966.4399, 966.4438, 966.4666, 
966.4139, 966.3791, 966.3496, 966.3203, 966.2857, 966.2732, 966.2787, 966.2605, 
966.2548, 966.2368, 966.2141, 966.2078, 966.1792, 966.1958, 966.1841, 966.1838, 
966.1778, 966.127, 966.1378, 966.1152, 966.1035, 966.126, 966.0864, 966.1371, 966.1143, 
966.0972, 966.1028, 966.1141, 966.1089, 966.1031, 966.1031, 966.1075, 966.051, 
966.1406, 966.1064, 966.1511, 966.1337, 966.1617, 966.1443, 966.1216, 966.0759, 
966.0251, 966.0699, 966.0417, 966.0471, 966.1147, 966.0694, 966.0919, 966.0938, 966.03, 
966.0525, 966.0794, 966.0342, 966.0557, 966.033, 966.014, 966.0384, 966.0267, 965.981, 
965.9422, 965.9099, 965.8817, 965.8833, 965.8833, 965.9128, 965.879, 965.8757, 
965.8965, 965.9114, 965.8699, 965.8426, 965.8389)

plot(data, type = "l")

find_peaks <- function (x, m){
  shape <- diff(sign(diff(x, na.pad = FALSE)))
  pks <- sapply(which(shape < 0), FUN = function(i){
  z <- i - m + 1
  z <- ifelse(z > 0, z, 1)
  w <- i + m + 1
  w <- ifelse(w < length(x), w, length(x))
  if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
  })
   pks <- unlist(pks)
   pks}

# 生成局部极值数据框（阈值为21分钟宽度）
maximums_df <- data.frame(trend = 'max', index = find_peaks(data, m = 10))
minimums_df <- data.frame(trend = 'min', index = find_peaks(-data, m = 10))
maxmin_df <- rbind(maximums_df, minimums_df)[order(rbind(maximums_df, 
             minimums_df)$index),]

# 多段趋势划分函数
segment_pressure_trend <- function(data, maxmin_df, time_res = 20, threshold = 0.1) {
  # 补充首尾端点
  full_points <- rbind(
    data.frame(trend = "start", index = 1),
    maxmin_df,
    data.frame(trend = "end", index = length(data))
  )
  
  segments <- list()
  
  # 遍历相邻点对
  for (i in 1:(nrow(full_points)-1)) {
    start_idx <- full_points$index[i]
    end_idx <- full_points$index[i+1]
    segment_data <- data[start_idx:end_idx]
    time_hours <- (end_idx - start_idx) * time_res / 60
    
    # 计算线性趋势斜率
    lm_fit <- lm(segment_data ~ seq_along(segment_data))
    slope <- coef(lm_fit)[2]
    rate_per_hour <- slope * (60 / time_res)  # 转换为毫巴/小时
    
    # 判断趋势类型
    if (abs(rate_per_hour) > threshold) {
      trend_type <- ifelse(rate_per_hour > 0, "上升", "下降")
    } else {
      trend_type <- "平稳"
    }
    
    segments[[i]] <- data.frame(
      start_idx = start_idx,
      end_idx = end_idx,
      trend = trend_type,
      rate_per_hour = round(rate_per_hour, 3)
    )
  }
  
  # 合并连续同趋势的区间
  merged_segments <- list(segments[[1]])
  for (seg in segments[-1]) {
    last_merged <- merged_segments[[length(merged_segments)]]
    if (seg$trend == last_merged$trend) {
      last_merged$end_idx <- seg$end_idx
      last_merged$rate_per_hour <- mean(c(last_merged$rate_per_hour, seg$rate_per_hour))
      merged_segments[[length(merged_segments)]] <- last_merged
    } else {
      merged_segments[[length(merged_segments)+1]] <- seg
    }
  }
  
  do.call(rbind, merged_segments)
}

# 执行划分并输出结果
result <- segment_pressure_trend(data, maxmin_df)
print(result)

# 可视化分段结果
plot(data, type = "l", main = "气压趋势分段")
for (i in 1:nrow(result)) {
  lines(x = result$start_idx[i]:result$end_idx[i], y = data[result$start_idx[i]:result$end_idx[i]], col = i+1, lwd = 2)
}
legend("topright", legend = paste(result$trend, "(", result$rate_per_hour, "mb/h)"), col = 2:(nrow(result)+1), lwd = 2)

代码说明 #

• segment_pressure_trend函数：核心处理逻辑，负责补充端点、计算区间趋势、合并同趋势段
• 时间分辨率参数time_res：对应数据的采样间隔（20分钟），用于转换趋势率单位
• 阈值参数threshold：自定义的趋势显著性阈值（±0.1毫巴/小时）
• 合并连续同趋势段：避免因局部小波动产生过多无效分段，确保结果符合实际观测

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