基于python和GIS实现Sen-mk长时间序列趋势性分析

import numpy as npimport matplotlib.pyplot as pltimport pymannkendall as mktime = np.arange(1993, 2023)  # 时间范围：1993-2022年# 模拟NDVI数据（含轻微上升趋势+随机噪声，贴近真实数据）np.random.seed(42)  # 固定随机种子，保证结果可复现trend = 0.02 * np.arange(len(time))  # 上升趋势（斜率0.02）noise = np.random.normal(0, 0.05, len(time))  # 随机噪声ndvi_data = 0.5 + trend + noise  # 使用pymannkendall库直接计算Sen斜率和MK检验结果# original_test方法返回结果包含：趋势方向、p值、Z值、Sen斜率等mk_result = mk.original_test(ndvi_data)# 提取关键结果trend_direction = mk_result.trend  # 趋势方向（'increasing'/'decreasing'/'no trend'）p_value = mk_result.p  # p值（与Z值等价，p<0.05表示显著）z_value = mk_result.z  # Z值（判断显著性）sen_slope = mk_result.slope  # Sen斜率（趋势速率）# ---------------------- 3. 结果解读 ----------------------print("=" * 50)print("SEN-MK长时间序列趋势分析结果")print("=" * 50)print(f"时间序列长度：{len(time)}个时间点")print(f"Sen斜率：{sen_slope:.4f}")print(f"趋势方向：{trend_direction}")print(f"MK检验Z值：{z_value:.4f}")print(f"MK检验p值：{p_value:.4f}")# 根据Z值判断显著性if abs(z_value) >= 2.58:    significance = "极显著（99%置信水平）"elif 1.96 <= abs(z_value) < 2.58:    significance = "显著（95%置信水平）"elif 1.645 <= abs(z_value) < 1.96:    significance = "弱显著（90%置信水平）"else:    significance = "不显著"print(f"趋势显著性：{significance}")print("=" * 50)# ---------------------- 4. 结果可视化（直观展示趋势）----------------------plt.rcParams['font.sans-serif'] = ['SimHei']  # 解决中文显示问题plt.figure(figsize=(12, 6))# 绘制原始时间序列plt.plot(time, ndvi_data, 'o-', color='#1f77b4', label='原始NDVI数据', linewidth=1.5, markersize=4)trend_line = sen_slope * (time - time[0]) + ndvi_data[0]plt.plot(time, trend_line, '-', color='#ff7f0e', linewidth=2, label=f'Sen趋势线（斜率={sen_slope:.4f}）')# 添加显著性标记plt.text(time[0] + 2, np.max(ndvi_data) - 0.05,          f'趋势：{trend_direction}\n显著性：{significance}\nSen斜率：{sen_slope:.4f}',         bbox=dict(boxstyle='round', facecolor='lightgray', alpha=0.8),         fontsize=10)# 图表美化plt.xlabel('年份', fontsize=12)plt.ylabel('NDVI植被指数', fontsize=12)plt.title('长时间序列NDVI数据SEN-MK趋势分析（1993-2022）', fontsize=14, pad=20)plt.legend(loc='best', fontsize=10)plt.grid(True, alpha=0.3)plt.tight_layout()# 保存图片（可选）plt.savefig('sen_mk_trend_analysis.png', dpi=300, bbox_inches='tight')plt.show()# ---------------------- 补充：手动实现Sen斜率（可选，理解原理）----------------------def calculate_sen_slope(data):    """手动计算Sen斜率"""    n = len(data)    slopes = []    # 遍历所有两两数据点，计算斜率    for i in range(n):        for j in range(i + 1, n):            slope = (data[j] - data[i]) / (j - i)            slopes.append(slope)    # 返回斜率的中位数    return np.median(slopes)# 手动计算Sen斜率（与pymannkendall结果一致）manual_sen_slope = calculate_sen_slope(ndvi_data)print(f"手动计算Sen斜率：{manual_sen_slope:.4f}")