NDVI时间序列逐像元计算Hurst指数

自相似性和长期依赖性是自然界普遍存在的现象，并在水文、气候、地质和地震等领域广泛运用，Hurst指数是描述该现象的有效方法。目前，Hurst指数的估算方法有多种，如绝对值法、聚合方差法、周期图法、小波分析法、残差分析法和R/S分析法等，有关研究表明：R/S分析法和小波分析法估算的Hurst指数要比其他方法估算的结果更可靠。基于重标极差(R/S)分析方法的Hurst指数最早是由英国水文学家Hurst在研究尼罗河水库流量和储存能力的关系时提出，本研究即采用R/S分析法。

NDVI 时间序列 NDVIi，i=1,2,3,4,…,n，对于任意正整数 m，定义该时间序列:对于比值R(m) / S(m)≌ R/S，若存在如下关系，则说明分析的时间序列存在Hurst现象，H称为Hurst指数。在双对数坐标系(lni,lnR/S) 中用最小二乘法拟合式得到。

R语言计算Hurst指数

R语言有专门针对水文学的一系列程序包（参考文献2），其中reservoir包中带有Hurst函数。使用terra包的栅格计算app函数，指定数据和Hurst函数，即可计算，在这里我使用了8核并行计算。具体代码如下：

library(terra)
library(reservoir)
#数据读取
GIMMS <- rast("./Yearmaxdata/YearMax.tif")  #读取GIMMS数据
#Hurst计算
GIMMS_hurst <- app(GIMMS, Hurst, cores=8)
#计算结果输出
writeRaster(GIMMS_hurst, filename = "./Yearmaxdata/GIMMS_hurst.tif")

Hurst指数计算结果

根据H的大小可以判断NDVI序列的持续性。Hurst指数取值包括三种形式：

• 若0.5<H<1，表明时间序列是一个持续性序列，具有长期相关的特征；
• 若H=0.5，则说明NDVI时间序列为随机序列，不存在长期相关性；
• 若0<H<0.5，表明NDVI时间序列具有反持续性。
• H值越接近于0，其反持续性越强；越接近1，其持续性越强。

下面是Hurst函数的R代码，供大家学习：

#' @title Hurst coefficient estimation
#' @description Hurst coefficient estimation.
#' @param Q   vector or annualized time series object. Net inflows or streamflow totals.
#' @return Returns an estimate of the Hurst coefficient, H (0.5 < H < 1).
#' @examples Q_annual <- aggregate(resX$Q_Mm3) #convert monthly to annual data
#' Hurst(Q_annual)
#' @import stats
#' @export
Hurst <- function(Q) {
    if (frequency(Q) != 1)
        warning("Q is not an annualized series")
    if (is.ts(Q) == FALSE && is.vector(Q) == FALSE)
        stop("Q must be time series or vector object")
    Q <- Q - mean(Q)
    max.Q <- max(cumsum(Q))
    min.Q <- min(cumsum(Q))
    RS <- (max.Q - min.Q) / sd(Q)
    H <- log(RS) / log(length(Q))
    return(H)
}

我试验了好几个Hurst有关的程序包，比如pracma和PerformanceAnalytics，计算一列数字可以，但是用于上面的栅格计算会报错。

如果觉得不错还请分享，在看

参考文献

1. https://www.jianshu.com/u/b3c93dcffd15
2. https://cran.r-project.org/web/views/Hydrology.html
3. https://github.com/cran/reservoir/blob/master/R/hurst_basic.R
4. 严恩萍, 林辉, 党永峰, 夏朝宗. 2000-2012年京津风沙源治理区植被覆盖时空演变特征[J]. 生态学报, 2014, 34(17): 5007-5020.