上一期讲了多元线性回归残差分析的计算，接下来讲解如何根据残差分析结果进行驱动因素的划分，计算贡献率等。
下表是驱动因素划分规则，根据上一期文章的残差分析结果，计算NDVI观测值、预测值、残差值的斜率，然后根据斜率划分驱动因素和计算贡献率。
| NDVI观测值slope | 重分类值 | NDVI预测值slope | 重分类值 | NDVI残差slope | 残差重分类值 | 驱动因素 | 驱动因素重分类值 |
|---|---|---|---|---|---|---|---|
| >0 | 1 | >0 | 1 | >0 | 3 | G&CC&HA | 3 |
| >0 | 1 | <0 | 4 | G&CC | 4 | ||
| <0 | 2 | >0 | 3 | G&HA | 6 | ||
| <0 | -1 | <0 | 2 | <0 | 4 | B&CC&HA | -8 |
| <0 | 2 | >0 | 3 | B&CC | -6 | ||
| >0 | 1 | <0 | 4 | B&HA | -4 |
驱动因素划分与面积占比计算
对于驱动因素的划分，实质上就是一个针对NDVI观测值、预测值、残差值斜率的重分类和乘法的组合。
NDVI观测值、预测值、残差值斜率计算
在这里先使用了Theil-Sen's median的方法计算了斜率，目的是得到NDVI观测值、预测值、残差值的趋势。注意，后来有公众号粉丝跟我反映，这个斜率计算使用一元线性回归更合适，建议改为一元线性回归计算斜率，可以参考再次条的文章，年底太忙了，没来得及改。
#计算NDVI观测值、预测值、残差值的趋势
library(terra)
#读取NDVI观测值
ndvi_obs = rast(fs::dir_ls("./GIMMSYM", regexp = ".tif"))
ndvi_obs
ndvi_pred = rast("./ResidualAnalysis/ndvi_pred.tif")
ndvi_pred
ndvi_res = rast("./ResidualAnalysis/ndvi_res.tif")
ndvi_res
#使用Sen+Mk方法计算斜率
#sen'slope计算公式
fun_sen = function(x){
if(length(na.omit(x)) <34) return(c(NA, NA, NA)) #删除数据不连续含有NA的像元
MK_estimate <- trend::sens.slope(ts(na.omit(x), start = 1982, end = 2015, frequency = 1), conf.level = 0.95) #Sen斜率估计
slope <- MK_estimate$estimate
MK_test <- MK_estimate$p.value
Zs <- MK_estimate$statistic
return(c(Zs, slope, MK_test))
}
#计算senslope
ndvi_obs_senslope = app(ndvi_obs, fun_sen, cores = 56)
names(ndvi_obs_senslope) = c("Z", "slope", "p-value")
ndvi_pred_senslope = app(ndvi_pred, fun_sen, cores = 56)
names(ndvi_pred_senslope) = c("Z", "slope", "p-value")
ndvi_res_senslope = app(ndvi_res, fun_sen, cores = 56)
names(ndvi_res_senslope) = c("Z", "slope", "p-value")
plot(ndvi_obs_senslope)
plot(ndvi_pred_senslope)
plot(ndvi_res_senslope)
植被变化驱动因素划分
按上面的表格，对计算得到的NDVI观测值斜率、NDVI预测值斜率、NDVI残差斜率进行重分类 对NDVI观测值、预测值、残差值斜率的重分类结果进行相乘，得到最后的驱动因子重分类值
#驱动因素判定及贡献率计算
#对Theil-Sen median趋势分析结果进行重分类，坡度>0增加赋值1，坡度<0减少赋值-1
slprcl = matrix(c(-20, 0, -1,
0, 20,1),
ncol = 3, byrow = T)
ndvi_obs_senslope_srcl = classify(ndvi_obs_senslope[["slope"]], slprcl, include.lowest = F, right = F)
#重分类NDVI预测值slope
predslope_rcl = matrix(c(-1, 0, 2,
0, 1,1),
ncol = 3, byrow = T)
ndvi_pred_senslope_rcl = classify(ndvi_pred_senslope[["slope"]], predslope_rcl, include.lowest = F, right = F)
plot(ndvi_pred_senslope_rcl)
#重分类NDVI残差slope
resslope_rcl = matrix(c(-1, 0, 4,
0, 1,3),
ncol = 3, byrow = T)
ndvi_res_senslope_rcl = classify(ndvi_res_senslope[["slope"]], resslope_rcl, include.lowest = F, right = F)
plot(ndvi_res_senslope_rcl)
ndvi_drivefactor = ndvi_obs_senslope_srcl*ndvi_pred_senslope_rcl*ndvi_res_senslope_rcl
plot(ndvi_drivefactor)
writeRaster(ndvi_drivefactor, filename = "./ResidualAnalysis/ndvi_drivefactor.tif")
将上面的结果在GIS软件中分类符号化，即可得到植被变化的空间分布图。具体可以查看#QGIS专题中的绘图教程。
驱动因素面积统计，本质上就是使用R计算了一下不同像元值的频率分布，下面给出代码：
## 统计不同驱动因素的面积占比
#读取气候驱动因子和对应汉语解释表格
library(tidyverse)
factortype = readxl::read_excel(path = "./ResidualAnalysis/factortype.xlsx", sheet = 1, col_names = T) %>%
mutate(drivefactor = factor(drivefactor, levels = c("-8", "-6", "-4", "-3", "3", "4", "6", "8")))
#统计不同驱动因素的占比
factorrast_df = values(ndvi_drivefactor) %>%
as.data.frame() %>%
set_names("drivefactor") %>%
na.omit() %>%
table() %>%
as.data.frame() %>%
mutate(Percent = Freq/sum(Freq)) %>%
left_join(factortype, by = "drivefactor")
factorrast_df
#输出计算结果
openxlsx::write.xlsx(factorrast_df, file = "./ResidualAnalysis/drivefactor.xlsx")
人类活动和气候变化贡献率计算
贡献率的计算非常简单，就是做一个斜率的除法，需要注意的是有些斜率相除的结果会出现负无穷等情况，需要使用terra包的clamp函数将贡献率的值域进行一个限定。
#气候变化和人类活动贡献率计算
#气候变化贡献率
ndvi_cc = clamp(ndvi_pred_senslope[["slope"]]/ndvi_obs_senslope[["slope"]], 0, 1, values = F) #clamp可以去除异常值
ndvi_cc
plot(ndvi_cc)
#人类活动贡献率
ndvi_ha = clamp(ndvi_res_senslope[["slope"]]/ndvi_obs_senslope[["slope"]], 0, 1, values = F) #clamp可以去除异常值
ndvi_ha
plot(ndvi_ha)
计算完成后，重复上面的绘图和像元频率计算步骤，即可得到贡献率的分布图和面积占比，在这不再重复，需要的同学可以自己绘制一下即可。
以上就完成了NDVI和气候因子（气温、降水）的多元线性回归残差分析的所有计算工作。觉得不错还请分享、点赞、再看三连。
点击阅读原文查看视频课程讲解，更多课程欢迎扫码加入培训班
