使用R语言将年日最大降雨量数据解聚为小时降雨数据的方法咨询

我们这里优先选用tempdisagg包完成操作，它的语法简单、对新手友好，无需额外的水文专业参数配置就能得到误差可控的解聚结果。如果你需要更贴合降雨水文特征的模拟结果，可以参考文末的BLRPM包实现示例。

前置准备：安装加载依赖包 #

# 安装包（首次运行执行即可）
install.packages("tempdisagg")
install.packages("tidyverse") # 用于数据处理和结果导出

# 加载包
library(tempdisagg)
library(tidyverse)

步骤1：导入你的降雨数据 #

你提供的原始数据可以直接用以下代码构造，如果你是从本地csv读取，替换为read.csv("你的文件路径.csv")即可：

rain_daily <- tibble(
  YEAR = c(1955,1956,1957,1958,1959,1960,1961,1962,1963,1964,1965,1966,1967,1968,1969,1970,1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004),
  MAX_DAILY = c(79.76,64.77,67.82,44.2,52.83,60.45,77.72,50.04,59.94,90.68,101.09,55.88,97.79,76.96,77.98,48.77,68,69.769,69.769,69.769,52.1,69.769,74.3,65.1,79.2,77.9,55.1,58.3,86.6,53.2,60.6,57.6,108.1,132.1,45.9,55.8,118.6,59.6,86,52,78.5,58,48.3,60.4,72.5,83.9,58.3,57.9,68.7,70.1)
)

步骤2：批量解聚所有年份的日降雨为小时数据 #

我们对每一年的日最大降雨量单独执行解聚，每个日值对应输出24个小时的降雨量：

# 设置随机种子保证结果可复现
set.seed(123)

# 批量处理所有年份
rain_hourly <- map_dfr(1:nrow(rain_daily), function(i){
  # 提取单年的日降雨值
  daily_val <- rain_daily$MAX_DAILY[i]
  year_val <- rain_daily$YEAR[i]
  
  # 构造低频时间序列（1个日值，解聚为24个小时值）
  ts_low <- ts(daily_val, frequency = 1)
  
  # 执行解聚，选用Denton方法，无偏且适合无辅助指标的解聚场景
  mod <- td(ts_low ~ 1, to = 24, method = "denton")
  hourly_val <- predict(mod)
  
  # 整理结果返回
  tibble(
    YEAR = year_val,
    HOUR = 1:24,
    HOURLY_RAIN = as.numeric(hourly_val)
  )
})

步骤3：结果校验与导出 #

# 校验：每一年24小时降雨总和应该等于原日最大降雨量，误差可忽略
check <- rain_hourly %>% 
  group_by(YEAR) %>% 
  summarise(SUM_HOURLY = sum(HOURLY_RAIN)) %>% 
  left_join(rain_daily, by = "YEAR")
head(check) # 查看前几行的校验结果

# 导出结果为csv文件
write.csv(rain_hourly, "小时降雨解聚结果.csv", row.names = F)

补充：BLRPM包（降雨专用解聚）实现示例 #

BLRPM是专门针对降雨时间序列解聚的包，模拟结果更贴合降雨的小时分布特征，基础实现代码如下：

# 安装加载包
install.packages("BLRPM")
library(BLRPM)

# 单年解聚示例，批量处理逻辑和上述tempdisagg的批量代码一致
# 以1955年79.76mm的日降雨为例
blrpm_mod <- BLRPM(D = 24, # 总历时24小时
                   Rtot = 79.76, # 总降雨量
                   # 以下为雨型参数，可根据你当地的小时降雨统计参数调整
                   a = 0.4, b = 1.2, c = 0.6)
hourly_rain_blrpm <- blrpm_mod$rain

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