基于R语言bfast的单变量森林覆盖时间序列趋势分析咨询

Hey there! 针对你手头的16年森林覆盖百分比时间序列数据，我整理了一套实用的R分析流程，帮你精准找出森林覆盖的增减时段，一步步来：

第一步：把数据整理成结构化格式 #

你给的原始数据是按年份分组罗列的，首先得转换成R能轻松处理的表格形式。我写了一段代码适配你的样例数据，后续你补全剩余年份的数据就行：

# 输入你的样例数据（把后面的...换成完整数据即可）
raw_data <- c(2000, rep(98,3), rep(99,9),
              2001, rep(98,9), rep(99,2),
              2002, 93,94, rep(95,8),96,
              2003, rep(96,5), rep(95,2), rep(96,4),
              2004, 55,54,55,55,53,49,46,45,52,59,65,
              2005,53,51,49,52,52,50,50,52,61,70,74,
              2006,32,33,46)

# 用tidyverse把数据转换成数据框
library(tidyverse)
# 先定位所有年份的位置
year_positions <- which(raw_data %in% 2000:2015) # 假设数据覆盖2000-2015共16年
data_frames <- list()

# 循环拆分每个年份的数据
for(i in seq_along(year_positions)){
  start_idx <- year_positions[i]
  end_idx <- if(i < length(year_positions)) year_positions[i+1]-1 else length(raw_data)
  current_year <- raw_data[start_idx]
  cover_values <- raw_data[(start_idx+1):end_idx]
  data_frames[[i]] <- data.frame(Year = current_year, Forest_Cover = cover_values)
}

# 合并成最终数据框
forest_df <- bind_rows(data_frames)
# 可选：转换成时间序列对象，方便后续分析
forest_ts <- ts(forest_df$Forest_Cover, start = min(forest_df$Year), 
                frequency = nrow(filter(forest_df, Year == 2000)))

第二步：先画个图直观看看趋势 #

先通过可视化快速把握整体变化，用年度均值展示趋势，同时保留原始数据点看离散情况：

ggplot(forest_df, aes(x = Year, y = Forest_Cover)) +
  geom_point(alpha = 0.6, color = "green") +
  geom_line(aes(group = Year), stat = "summary", fun = "mean", 
            color = "darkgreen", size = 1.2) +
  labs(title = "森林覆盖百分比年度变化趋势", 
       x = "年份", y = "森林覆盖百分比(%)") +
  theme_minimal()

第三步：量化分析趋势与识别增减时段 #

3.1 线性回归看整体趋势 #

先跑个简单线性回归，快速判断整体是增长还是下降：

lm_fit <- lm(Forest_Cover ~ Year, data = forest_df)
summary(lm_fit)
# 看输出里Year的系数：正数代表整体增长，负数代表整体下降，p值判断显著性

3.2 分段回归找趋势突变点 #

这是核心步骤——用分段回归找出森林覆盖趋势发生变化的年份，明确不同的增减时段。用segmented包实现：

install.packages("segmented")
library(segmented)

# 基于线性回归模型构建分段回归
seg_fit <- segmented(lm_fit, seg.Z = ~Year)
summary(seg_fit)

# 可视化分段回归结果，直观看到突变点和不同时段的趋势
plot(seg_fit, res = TRUE, main = "森林覆盖分段回归趋势", 
     xlab = "年份", ylab = "森林覆盖百分比(%)")

输出里的Breakpoint就是趋势发生突变的年份，比如从持续增长突然转为下降的节点。

3.3 Mann-Kendall非参数检验（稳健备选） #

如果你的数据不满足正态分布，这个非参数方法更靠谱，能稳健检验趋势的显著性：

install.packages("Kendall")
library(Kendall)

mk_result <- MannKendall(forest_ts)
print(mk_result)
# tau值为正代表上升趋势，负代表下降；p值<0.05说明趋势显著

第四步：划分时段并总结变化 #

根据分段回归得到的突变点，把16年划分成不同时段，然后计算每个时段的趋势斜率，明确是增长还是下降：

# 提取突变点（取整数年份）
breakpoints <- floor(seg_fit$psi[, "Est."])
# 构建时段区间
periods <- data.frame(
  起始年份 = c(min(forest_df$Year), breakpoints),
  结束年份 = c(breakpoints - 1, max(forest_df$Year))
)

# 计算每个时段的趋势
period_trends <- periods %>%
  rowwise() %>%
  mutate(
    时段数据 = list(filter(forest_df, Year >= 起始年份 & Year <= 结束年份)),
    趋势模型 = list(lm(Forest_Cover ~ Year, data = 时段数据)),
    斜率 = coef(趋势模型)[["Year"]],
    变化方向 = ifelse(slope > 0, "增长", "下降")
  ) %>%
  select(起始年份, 结束年份, 斜率, 变化方向)

print(period_trends)

运行这段代码后，你就能清晰看到每个时段森林覆盖的变化方向和幅度啦。

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