如何用ggpmisc或其他R包识别水文时间序列的真实峰值？

我正在处理一段为期两周、间隔15分钟的河流量观测时间序列，想要确定该时段内河流量峰值对应的dateTime。尝试过ggpmisc::stat_peaks()、pracma::findpeaks()和cardidates::peakwindow()等工具，但均未得到符合人工目视判断的真实峰值结果。以下是问题复现及解决方案：

数据获取与预处理 #

library(tidyverse)
library(dataRetrieval)
library(ggpmisc)

# 下载河流量时间序列数据
dat <- readNWISuv(siteNumbers = "01576381",
                  parameterCd = "00060",
                  startDate = "2024-04-01",
                  endDate = "2024-04-15",
                  tz = "America/New_York") %>% 
  renameNWISColumns() %>% 
  rename(flow = Flow_Inst) %>% 
  select(dateTime, flow)

# 查看数据结构
head(dat)
#>              dateTime flow
#> 1 2024-04-01 00:00:00 25.0
#> 2 2024-04-01 00:15:00 25.6
#> 3 2024-04-01 00:30:00 25.0
#> 4 2024-04-01 00:45:00 25.6
#> 5 2024-04-01 01:00:00 25.1
#> 6 2024-04-01 01:15:00 25.6

原始时间序列可视化 #

dat %>% 
  ggplot(aes(x = dateTime, y = flow)) +
  geom_line() +
  theme_bw()

stat_peaks的尝试结果 #

默认参数识别 #

dat %>% 
  ggplot(aes(x = dateTime, y = flow)) +
  geom_line() +
  stat_peaks(color = 'red') +
  theme_bw()

开启strict参数识别 #

dat %>% 
  ggplot(aes(x = dateTime, y = flow)) +
  geom_line() +
  stat_peaks(color = 'red',
             strict = TRUE) +
  theme_bw()

理想识别结果（人工目视判断） #

解决方案 #

1. 优化ggpmisc::stat_peaks参数 #

默认参数会识别所有局部极小值，导致大量高频小峰值。通过设置窗口范围与高度阈值，过滤掉无意义的微小波动：

可视化层面调整

dat %>% 
  ggplot(aes(x = dateTime, y = flow)) +
  geom_line() +
  stat_peaks(color = 'red',
             span = 96, # 窗口半宽，对应24小时（96个15分钟间隔点）
             ignore_threshold = 0.5, # 峰值与窗口内最小值的最小高度差
             strict = TRUE) +
  theme_bw()

• span：定义判断峰值的滑动窗口半宽，单位为数据点数量。设置为96可过滤日内小波动，聚焦跨日的真实峰值。
• ignore_threshold：设定峰值必须高于窗口内最小值的阈值，过滤微小起伏。

提取峰值数据

若需直接获取峰值对应的dateTime与流量值，使用底层的ggpmisc::find_peaks函数：

peaks_dat <- dat %>% 
  mutate(is_peak = ggpmisc::find_peaks(flow, span = 96, ignore_threshold = 0.5, strict = TRUE)) %>% 
  filter(is_peak)

print(peaks_dat)

2. 推荐其他工具方案 #

pracma::findpeaks

该函数支持更灵活的峰值规则定义，比如最小峰值高度、峰值间最小间隔：

library(pracma)

# 提取流量向量
flow_vals <- dat$flow

# 识别峰值：设置最小峰值高度、峰值间最小间隔（96个点=24小时）
peak_info <- findpeaks(flow_vals, minpeakheight = 26, minpeakdistance = 96)

# 提取对应数据
peaks_pracma <- dat[peak_info[, 2], ]

forecast::findPeaks

用法简洁，通过阈值与窗口过滤小峰值：

library(forecast)

# 识别峰值
peak_indices <- findPeaks(flow_vals, threshold = 0.5, span = 96)

# 提取对应数据
peaks_forecast <- dat[peak_indices, ]

核心思路 #

所有工具的核心逻辑一致：通过定义窗口范围过滤高频波动，通过高度阈值区分真实峰值与噪声。需根据数据的波动特征调整参数，比如窗口大小对应你认为的“真实峰值”时间间隔，高度阈值对应峰值与背景波动的最小差值。

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