Hurst指数
Hurst指数（Hurst Exponent）是衡量时间序列长期记忆性和趋势持续性的重要统计工具，由英国水文学家Harold Edwin Hurst于1951年提出。H值的范围通常在0到1之间，其中H=0.5表示随机游走（即无记忆性），H>0.5表示存在持久性趋势，而H<0.5则表明存在反持久性或均值回复特征。
Hurst指数主流计算方法：
- 重标极差分析法（R/S Analysis） 步骤 ① 将序列分割为子区间 ② 计算每个子区间的极差（Range）与标准差（Standard Deviation） ③ 通过线性回归拟合 log(R/S)与log(n)的斜率即为H值 优势：对非高斯分布数据稳健 局限：对短序列敏感，易高估H值
- 去趋势波动分析（DFA, Detrended Fluctuation Analysis） 通过消除局部趋势干扰，更适合含非平稳噪声的序列（如生理信号、金融数据）。
- 小波变换法 利用多尺度分析捕捉序列的局部特征，适用于高频数据。
Hurst指数的主要应用场景：
金融市场分析：Hurst指数广泛应用于金融市场的波动性研究中，帮助投资者理解市场趋势是随机、持续还是具有回归均值的倾向。可以帮助分析股票的买卖点。
网络流量预测：在网络工程领域，H指数可用于分析和预测网络流量模式，有助于规划资源分配和提高服务质量。
气候科学：对于气候变化研究，Hurst指数可以帮助科学家们识别出温度、降水等气象要素的时间序列中的长期趋势及其变化规律。
地质与环境研究：例如，在地下水位变动、河流流速等自然过程的研究中，H指数能够提供对这些过程长期行为的理解。
代码实现
我这里有2000年至2023年NDVI指数的时间序列数据。计算时间序列栅格数据Hurst指数其实就是计算每个像元的Hurst指数，最终得到Hurst指数的栅格文件。
在我之前的博客中利用了python进行Theil-Sen Median斜率估计和Mann-Kendall趋势分析。这里计算Hurst指数的实现思路相同，所以我在之前的代码基础上进行了修改，使得可以同时分析Hurst指数（持续性判断）、Sen趋势（趋势方向判断）、MK显著性（显著性判断）。
借助Deepseek-R1模型帮我完善了代码，完整代码如下：
代码逻辑简要分析：calculate_hurst函数就是计算Hurst指数的实现，采用了重标极差分析法（R/S Analysis）。时间序列长度大于10的才会分析，当然数据量越大越好。classify_trend函数根据不同的H值、Slope值和Z值划分类别。_encode_class这是一个私有函数，对划分的类别进行编码，不同的数字代码不同的情况，如数字1表示极显著持续增加。process_pixel函数就是处理单个像元，包括计算Hurst指数、Theil-Sen斜率及MK显著性。main函数首先是匹配所有的tif文件并按照年份排序，然后构建像元数据序列并进行并行处理，最后将结果输入到新的tif文件，这个tif文件包含四个波段，分别是编码波段（即1-10）、Theil-Sen斜率波段、Hurst指数波段、MK显著性波段。
import glob
import re
import numpy as np
import rasterio
from pymannkendall import original_test
from scipy.stats import theilslopes
from numba import jit
from tqdm import tqdm
# ======================
# 核心计算函数
# ======================
@jit(nopython=True, cache=True)
def calculate_hurst(ts):
"""优化Hurst指数计算（R/S方法）"""
n = len(ts)
if n < 10 or np.all(np.isnan(ts)):
return np.nan
# 数据清洗
ts = ts[~np.isnan(ts)]
n = len(ts)
if n < 10:
return np.nan
# 累积离差序列
mean_ts = np.mean(ts)
deviations = ts - mean_ts
cumulative_dev = np.cumsum(deviations)
# 计算极差
R = np.max(cumulative_dev) - np.min(cumulative_dev)
S = np.std(ts)
if S == 0 or R == 0:
return np.nan
return np.log(R / S) / np.log(n)
def classify_trend(slope, z, H):
"""直观趋势分类系统"""
if np.isnan(slope) or np.isnan(z) or np.isnan(H):
return 0 # 无效数据
# 趋势方向判断
if slope > 0.001:
trend = 1 # 增加
elif slope < -0.001:
trend = -1 # 减少
else:
return 10 # 无显著变化
# 显著性判断
if abs(z) > 2.58:
sig = 3 # 极显著
elif abs(z) > 1.96:
sig = 2 # 显著
elif abs(z) > 1.65:
sig = 1 # 微显著
else:
return 10 # 无显著变化
# 持续性判断
if H > 0.65:
pers = 3 # 强持续
elif H > 0.55:
pers = 2 # 可能持续
elif H < 0.45:
pers = -2 # 可能反转
else:
pers = 1 # 随机波动
# 生成组合编码
return _encode_class(trend, sig, pers)
def _encode_class(trend, sig, pers):
"""编码逻辑（私有函数）"""
if trend == 1:
if sig == 3:
if pers == 3: return 1
if pers == 2: return 2
if pers == -2: return 8
elif sig == 2:
if pers == 3: return 3
elif sig == 1:
if pers == 1: return 6
elif trend == -1:
if sig == 3:
if pers == 3: return 4
if pers == 2: return 9
if pers == -2: return 5
elif sig == 2:
if pers == 1: return 7
return 99 # 未定义类别
# ======================
# 数据处理流程
# ======================
def process_pixel(ts):
"""处理单个像素时间序列"""
ts = np.array(ts, dtype=np.float32)
valid_ts = ts[~np.isnan(ts)]
# 有效性检查
if len(valid_ts) < 10:
return (0, np.nan, np.nan, np.nan)
try:
# Theil-Sen斜率
slope, *_ = theilslopes(valid_ts, np.arange(len(valid_ts)))
# Mann-Kendall检验
mk_result = original_test(valid_ts)
# Hurst指数
H = calculate_hurst(valid_ts)
# 分类编码
code = classify_trend(slope, mk_result.z, H)
return (code, slope, H, mk_result.z)
except:
return (0, np.nan, np.nan, np.nan)
# ======================
# 主程序
# ======================
def main():
# 加载输入文件
tif_files = glob.glob('./NDVI_timeSeries/*.tif')
sorted_files = sorted(tif_files,
key=lambda x: int(re.search(r'\d{4}', x).group()))
# 初始化数据立方体
with rasterio.open(sorted_files[0]) as src:
meta = src.meta.copy() # 确保复制元数据
meta.update(count=4, dtype=rasterio.float32)
height, width = src.shape
num_pixels = height * width
# 构建三维数据阵列
data_cube = []
for f in tqdm(sorted_files, desc="加载数据"):
with rasterio.open(f) as src:
data_cube.append(src.read(1).astype(np.float32).flatten())
data_cube = np.array(data_cube).T # 转换为（像素数 × 时间步长）
# 并行处理像素
results = []
for pixel_ts in tqdm(data_cube, desc="处理像素", total=num_pixels):
results.append(process_pixel(pixel_ts))
# 重构结果矩阵
codes = np.array([r[0] for r in results]).reshape(height, width)
slopes = np.array([r[1] for r in results]).reshape(height, width)
hurst = np.array([r[2] for r in results]).reshape(height, width)
z_scores = np.array([r[3] for r in results]).reshape(height, width)
# 写入输出文件
with rasterio.open('Trend_Analysis_Result.tif', 'w', **meta) as dst:
dst.write(codes.astype(np.float32), 1)
dst.write(slopes.astype(np.float32), 2)
dst.write(hurst.astype(np.float32), 3)
dst.write(z_scores.astype(np.float32), 4)
dst.update_tags(
BAND_DESCRIPTIONS="Trend_Code, Slope, Hurst, Z_Score",
CLASS_LEGEND="""
1: 极显著持续增加 4: 极显著持续减少
2: 极显著可能持续增加 5: 极显著可能反转减少
3: 显著持续增加 9: 极显著可能持续减少
6: 微显著随机波动增加 7: 显著随机波动减少
8: 极显著可能反转增加 10: 无显著变化 0: 无效数据"""
)
if __name__ == "__main__":
main()
"""
1 : 极显著持续增加（H>0.65, Z>2.58, Slope>0.001）
2 : 极显著可能持续增加（H>0.55, Z>2.58, Slope>0.001）
4 : 极显著持续减少（H>0.65, Z<-2.58, Slope<-0.001）
5 : 极显著可能反转减少（H<0.45, Z<-2.58, Slope<-0.001）
3 : 显著持续增加（H>0.65, 1.96<Z<2.58, Slope>0.001）
9 : 极显著可能持续减少（H>0.55, Z<-2.58, Slope<-0.001）
6 : 微显著随机波动增加（0.45<H<0.55, 1.65<Z<1.96, Slope>0.001）
7 : 显著随机波动减少（0.45<H<0.55, Z<-1.96, Slope<-0.001）
8 : 极显著可能反转增加（H<0.45, Z>2.58, Slope>0.001）
10: 无显著变化
0 : 无效数据/数据不足
"""
# import numpy as np
# import matplotlib.pyplot as plt
# from matplotlib.colors import ListedColormap, BoundaryNorm
# import rasterio
#
# # 创建专属11色色阶
# cmap_11 = ListedColormap(
# plt.get_cmap('tab20')(np.linspace(0, 1, 11)), # 从tab20中抽取11个颜色
# name='Custom11'
# )
#
# # 定义数值-颜色映射规则
# class_boundaries = np.arange(-0.5, 11, 1) # 生成 [-0.5, 0.5, ..., 10.5]
# norm = BoundaryNorm(class_boundaries, cmap_11.N)
#
# with rasterio.open('Trend_Analysis_Result.tif') as src:
# codes = src.read(1)
#
# # 应用新色阶和映射规则
# plt.imshow(codes, cmap=cmap_11, norm=norm)
#
# # 配置精准colorbar
# cbar = plt.colorbar(
# ticks=np.arange(0, 11), # 显示0-10的整数刻度
# spacing='uniform', # 色块等间距
# boundaries=np.arange(-0.5, 11) # 严格对齐颜色边界
# )
# cbar.set_label('Trend Class', fontsize=12)
# cbar.ax.set_yticklabels([f'class {str(i)}' for i in range(11)]) # 强制显示整数标签
#
# plt.title('Vegetation Trend Analysis (2000-2023)', fontweight='bold')
# plt.show()
结果展示
上述代码最后用matpoltlib包（主要是需要调节一下colorbar颜色与类别相对应，需要11种颜色）绘制了编码波段的情况如下：
可以看到大部分都是类别1（极显著持续增加（H>0.65, Z>2.58, Slope>0.001））和类别10（无显著变化），还有少量的类别4（极显著持续减少）。
当然也可以从上述四个波段的tif文件中提取单个波段生成单波段的tif文件，代码如下：
这样就只有编码波段了，方便后续展示与处理。
import rasterio
def extract_single_band(input_path, output_path, band_number=1):
"""
从多波段TIFF文件中提取指定波段
参数：
input_path : str - 输入文件路径
output_path : str - 输出文件路径
band_number : int - 需要提取的波段号（从1开始）
"""
with rasterio.open(input_path) as src:
# 校验波段数量
if src.count < band_number:
raise ValueError(f"输入文件只有{src.count} 个波段，无法提取第{band_number}波段")
# 读取目标波段数据
band_data = src.read(band_number)
# 构建元数据
meta = src.meta.copy()
meta.update({
'count': 1,
'dtype': band_data.dtype,
'nodata': src.nodatavals[band_number - 1] if src.nodatavals else 0
})
# 写入输出文件
with rasterio.open(output_path, 'w', **meta) as dst:
dst.write(band_data, 1)
# 保留原始地理信息
dst.transform = src.transform
dst.crs = src.crs
# 添加波段描述
if band_number == 1:
dst.set_band_description(1, "Trend_Classification")
dst.update_tags(**src.tags())
# 使用示例
if __name__ == "__main__":
input_file = "Trend_Analysis_Result.tif" # 四波段文件路径
output_file = "Trend_Class.tif" # 输出文件路径
try:
extract_single_band(input_file, output_file)
print(f"成功提取趋势分类波段到：{output_file}")
except Exception as e:
print(f"处理失败：{str(e)}")
扫一扫
3万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
