R语言实战：用Terra包搞定NDVI的Sen+MK趋势分析（附完整代码）

最新推荐文章于 2026-03-12 00:37:58 发布

原创

最新推荐文章于 2026-03-12 00:37:58 发布 · 1.1k 阅读

R语言实战：用Terra包搞定NDVI的Sen+MK趋势分析（附完整代码）

如果你正在处理长时间序列的NDVI数据，想弄清楚植被覆盖到底是在变好还是变差，而且希望这个过程既快又准，那你来对地方了。Sen斜率估计结合Mann-Kendall检验（简称Sen+MK）是生态、水文和气候研究中的经典方法，它能稳健地评估趋势并检验其显著性。但面对动辄几十年的全国乃至全球栅格数据，逐像元计算简直就是一场噩梦，跑个几天几夜是常事。

今天，我们不谈复杂的理论推导，直接上手实战。我会带你用R语言中的terra包，配合trend包，构建一个高效、可复现、且支持并行计算的完整分析流程。从数据准备、并行计算函数编写，到结果可视化和制图，每一步都有清晰的代码和解释。无论你是研究生正在处理毕业论文数据，还是研究人员需要分析区域生态变化，这套方法都能让你事半功倍。

1. 环境准备与数据理解

在开始写代码之前，确保你的R环境已经就绪。我们主要依赖两个核心包：terra用于处理栅格数据，它比老的raster包速度更快，内存效率更高；trend包则提供了sens.slope函数，专门用于计算Sen斜率和MK检验。

# 安装必要的包（如果尚未安装）
install.packages("terra")
install.packages("trend")

# 加载包
library(terra)
library(trend)

你的数据应该是一系列时间序列的栅格文件，例如每年一幅的NDVI年均值图。文件命名最好有规律，比如NDVI_2000.tif， NDVI_2001.tif... 并放在同一个文件夹里。数据质量至关重要，要提前检查是否存在大量NA值（如云覆盖、水体），这会影响趋势计算的可靠性。

注意：terra包处理大型栅格数据集优势明显，但其语法与raster包略有不同。如果你熟悉raster，需要一点时间来适应，但转换后的效率提升是值得的。

2. 构建核心计算函数

Sen+MK分析的本质是对每个像元的时间序列值进行趋势计算。我们需要自己编写一个函数，封装sens.slope的计算，并处理好数据不完整的情况。

# 定义Sen+MK计算函数
fun_sen_mk <- function(x, start_year = 1982, end_year = 2015) {
  # x: 一个像元的时间序列数值向量
  # start_year, end_year: 时间序列的起始和结束年份
  
  # 1. 数据有效性检查
  # 剔除NA值，并检查剩余有效数据的长度
  x_clean <- na.omit(x)
  n_valid <- length(x_clean)
  
  # 如果有效数据太少（例如少于总年份的2/3），则返回NA
  total_years <- end_year - start_year + 1
  if (n_valid < 0.67 * total_years) {
    return(c(Z = NA, slope = NA, p_value = NA))
  }
  
  # 2. 创建时间序列对象并计算
  # 使用trend包的sens.slope函数
  ts_data <- ts(x_clean, start = start_year, end = end_year, fre

标签