EEMD原理介绍
集合经验模态分解（Ensemble Empirical Mode Decomposition，EEMD）可以用于研究NDVI和气候要素的非线性变化特征，得到不同时间尺度的模态函数及能反映时间序列数据长期变化趋势的残差分量。
EEMD方法通过只利用局部极值信息的"筛选"过程，将非线性和非平稳的时间序列数据分解为几个频率递减的内在模式函数（Intrinsic Mode Function, IMF）（）和长期趋势（），为了缓解"模态混叠"问题，提高分解结果的稳健性，EEMD在原始数据中加入一定数量的白噪声，然后取集合平均。集合经验模态分解的步骤包括：
第一步，添加一定振幅的高斯白噪声到原始时间序列数据
第二步，利用三次样条函数分别连接所有局部最大值和最小值，形成时间序列数据的上包络线和下包络线，然后利用时间序列数据减去上下包络线的均值。
第三步，确定是否满足在任何地方足够接近于零，满足则停止筛选，否则作为新的时间序列重复第二步，直到第k个包络的平均值满足标准，这样成为第一个。
第四步，从中减去得到余数，若余数包含振荡分量，则余数作为新的时间序列数据重复第二和第三步，直到余数变为单调函数或最多包含一个极值，最终时间序列被分解为数量为n的IMF分量和一个残差趋势分量。
第五步，重复l次第一至四步，每次添加相同振幅的高斯白噪声序列，最终将集合IMFs分量和趋势分量的平均值作为最终结果。在本文中高斯白噪声数量设置为100次，噪声振幅为原始时间序列数据的0.2倍标准差。
不同IMF分量的周期可以使用分量的长度和极值数量确定，公式如下：
式中，为IMF分量长度，为第个IMF分量极值的数量。
利用方差贡献率，可以评估不同时间尺度IMF分量的相对重要性，方差贡献率数值越大，说明对应的IMF分量越重要：
式中，，分别为第n个IMF分量的方差贡献率和方差，，分别对应为残差趋势分量的方差贡献率和方差。算法使用R语言Rlibeemd包、signal包编程实现。
EEMD编程实现
接下来讲一下如何把上面的原理转换为实际的代码实现。
数据预处理
实验数据来源于【数据共享】中国8km分辨率逐月NDVI数据集（1982-2022年）
library(terra)
library(tidyverse)
#读取逐月GIMMS3g+ NDVI月最大值合成数据
g3pmvclist = fs::dir_ls(path = "./g3pmvc/", regexp = ".tif$")
g3pmvc = rast(g3pmvclist)
#月最大值合成年度均值
yearIndex = sort(rep(1:41, 12))
yearIndex
g3pym = tapp(g3pmvc, yearIndex, mean, cores = 12)
names(g3pym) = seq(1982, 2022, 1)
g3pym
plot(g3pym)
把上面的栅格时间序列，使用全局统计的方法，转换为数据框时间序列：
#将栅格转为全局时间序列
g3pym_ts = global(g3pym, mean, na.rm = T, cores = 12)
EEMD计算与绘图
下面的代码完成了EEMD的计算与绘图，得到集合经验模态分解结果。
下面的代码把EEMD计算和绘图的过程打包成了一个函数，只要把计算的时间序列数据框输入到函数中，即可完成EEMD的计算和绘图工作，非常方便。
#EEMD计算
library(Rlibeemd)
library(signal)
library(ggpmisc)
#将EEMD计算，分量周期，方差贡献率，绘图组成一个大的函数
fun_eemd = function(input_ts, input_var){
#EEMD趋势分析
library(Rlibeemd)
input_ts_eemd = eemd(input_ts$mean, ensemble_size = 100, noise_strength = 0.2)
plot(input_ts_eemd)
input_ts_eemd = ts(input_ts_eemd, start = 1982, end = 2022, frequency = 1)
library(signal)
# 假设 emd_result 包含经验模态分解的结果，是一个矩阵，每一行是一个 IMF
# 分析每个 IMF 的频谱
for (i in 1:4) {
imf <- input_ts_eemd[, i]
# 计算频谱
spec <- spec.pgram(imf)
# 寻找主要频率或频谱峰值等周期信息
peak_freq <- spec$freq[which.max(spec$spec)]
#计算方差贡献率
imf_vars = diag(var(input_ts_eemd[,1:4])) #计算IMF分量的方差
res_var = stats::var(input_ts_eemd[,5]) #计算残差方差
#计算方差贡献率
total_var = sum(imf_vars)+res_var
imf_con = imf_vars/total_var #IMF分量方差贡献率
res_con = res_var/total_var #残差趋势分量的方差贡献率
# 输出结果
cat("IMF", i, "的主要频率（周期）为:", 1 / peak_freq, "年，方差贡献率为", imf_con[i]*100, "%。\n")
}
print(paste("残差趋势分量方差贡献率为", res_con*100, "%"))
eemd_result = input_ts_eemd %>%
as.data.frame() %>%
mutate(Year = seq(1982,2022,1)) %>%
pivot_longer(-Year, names_to = "type", values_to = "value") %>%
bind_rows(
input_ts %>%
mutate(Year = seq(1982,2022,1)) %>%
mutate(type = input_var) %>%
set_names(c("value", "Year", "type"))
) %>%
mutate(type = factor(type, levels = c(input_var, "IMF 1", "IMF 2", "IMF 3", "IMF 4", "Residual")))
eemd_result
# 将标签列表传递给facet_wrap中的labeller参数
p20 <- ggplot(eemd_result, aes(x = Year, y = value)) +
geom_line() +
facet_wrap(~type, ncol = 2, scales = "free_y") +
labs(x = "Year", y = "")+
theme_bw()+
theme(text = element_text(family = "serif"))
# 过滤出input_var对应的数据，且没有缺失值
input_data <- eemd_result %>%
dplyr::filter(type == input_var) %>% #filter这个函数在好几个包中出现，必须指定是dplyr中的，否则报错
drop_na(value)
# 在第一个分面中添加回归方程标注
p20 = p20 +
geom_smooth(data = input_data, method = "lm", formula = y ~ x, linewidth = 0.5)+
stat_poly_eq(data = input_data, aes(label = paste(after_stat(eq.label), after_stat(rr.label), after_stat(p.value.label), sep = "*\", \"*")),
formula = y ~ x, parse = TRUE, family = "serif", size = 2)
# 显示图形
p20
ggsave(filename = paste0("./eemd/", deparse(substitute(input_ts)),".jpg"), p20, width = 14, height = 8, dpi = 600, units = "cm")
}
fun_eemd(g3pym_ts, "NDVI")
参考文献
孙锐. 中国北方地区植被指数时空变化及其影响因素研究[D/OL]. 中国科学院大学（中国科学院地理科学与资源研究所）, 2021. http://ir.igsnrr.ac.cn/handle/311030/166303. https://github.com/helske/Rlibeemd 【数据共享】中国8km分辨率逐月NDVI数据集（1982-2022年） STOTEN插图复现|如何使用loess方法拟合NDVI时间序列变化
