目录

1 Mann-Kendall趋势检验算法

1.1 基本概念

1.2 具体原理

1.3 计算步骤

2 Python环境配置

3 Python完整代码

4 Python代码运行结果

---

1 Mann-Kendall趋势检验算法

1.1 基本概念

Mann-Kendall检验是一种非参数检验（无分布检验），用于检测时间序列数据中的单调趋势，不受数据分布限制，其优点是不要求样本遵从一定的分布，也不受少数异常值的干扰。该方法由Mann(1945)提出，后经Kendall(1975)完善，基于秩次比较原理，通过计算数据序列中所有可能配对的差值符号，得出趋势统计量及其显著性水平。常用于对降水、径流、气温和水质等要素时间序列变化趋势和突变点分析。

Mann-Kendall非参数检验方法常用于水质、径流量、温度、降水等水文气象时间序列变化趋势的显著性检验。该方法不需要测量值服从正太分布，不受缺失值和异常值的影响。在遥感中，MK趋势检验常用于分析植被指数(如NDVI、EVI)的长时间序列变化，帮助研究者识别植被覆盖变化的方向、速率及显著性水平，为生态环境监测、气候变化研究提供科学依据。

趋势分析可用于各种应用，例如水文学、气候学、城市测绘、农作物监测等，Mann-Kendall (MK) 趋势检验和 Sen 的斜率估计等非参数方法能够提供时间和空间尺度数据的趋势及其显著性水平。

Mann-Kendall检验的假设：原假设——时间序列中不存在趋势变化，备择假设——时间序列中存在趋势变化。其基本思想是通过比较每对数据点之间的大小关系来检查序列中的趋势，然后根据秩和的正负性来确定趋势的方向。

1.2 具体原理

对任意待检序列Xt(t=1,2,…,n)，n为待检序列长度，可定义统计量S：

其中：Xj和Xk为时间序列相应年份数据；n为时间序列长度；sgn(Xj-Xk)为符号函数。

当Xj＞Xk时，sgn值为1；当Xj＜Xk时，sgn值为-1；当Xj＝Xk时，sgn值为0。 

当n≥10时，统计量S近似服从正态分布，其期望和方差分别为： 

按照下式可构造标准化的检验统计量Z：

Z服从标准正态分布，当Z>0时，存在上升的趋势；当Z<0时，存在下降的趋势。对于给定的显著性水平α，如果|Z|≥Z1-α/2，说明时间序列存在显著向上或向下的趋势。

当|Z|≥1.64、1.96、2.58时则说明该时间序列分别通过了置信水平90%、95%、99%的显著性检验。

1.3 计算步骤

Mann-Kendall检验的计算过程包括以下步骤：

（1）对于长度为n的时间序列X，计算n(n-1)/2个变量关系差（d）:

（2）将所有d[i,j]的符号相乘，得到Sgn。即：

（3）计算所有d[i,j]的绝对值的秩和，即：

（4）计算标准正态分布的z统计量，即：

（5）根据显著性水平选择拒绝原假设的临界值，如果z值小于临界值，则接受原假设，否则拒绝原假设。

在Mann-Kendall检验确定存在趋势后，可以使用Sen's slope来计算趋势的斜率。Sen's slope是一种计算趋势斜率的方法，它利用中位数差来估计趋势的斜率，与线性回归方法不同，Sen's slope不受异常值的影响，因此更适用于含有异常值的数据。

2 Python环境配置

本代码使用的Python库包括：在基础数据处理库中，使用os库进行文件和目录操作，glob库用于批量查找和匹配TIF文件；在科学计算和数据分析库中，使用NumPy库进行数组操作和处理栅格数据，使用pymannkendall库实现Mann-Kendall趋势检验和Sens斜率计算；在栅格数据处理方面，使用rasterio库进行栅格数据的读写操作和元数据处理；在数据可视化方面，使用matplotlib库绘制累积距平图和MK统计量图。

# 基础数据处理
import os  # 文件路径管理和文件操作
import glob  # 批量匹配和查找TIF文件
import numpy as np  # 科学计算，数组操作和栅格数据处理

# 趋势检验与分析
import pymannkendall as mk  # 提供Mann-Kendall趋势检验和Sens斜率计算方法

# 栅格数据处理
import rasterio  # 读写栅格文件和处理元数据信息

# 数据可视化
import matplotlib.pyplot as plt  # 绘制累积距平图和MK统计量图
from matplotlib import rcParams  # 设置中文字体和图表样式

3 Python完整代码

代码主要功能步骤包括：首先读取NDVI栅格数据文件并按时间顺序排序；接着创建三维数组存储所有时间序列数据；然后定义趋势分类函数，基于Mann-Kendall检验的z值和slope值将趋势分为显著增加、微显著增加、稳定不变、微显著下降和显著下降五类；随后对每个像素点逐一进行MK检验并分类；同时计算全局平均NDVI的时间序列及其累积距平值和逐步MK检验统计量；绘制累积距平图和MK统计量图以直观展示变化趋势；最后将逐个像元的趋势检验结果保存为GeoTIFF格式栅格文件。

# ============================================================================
# 代码名称：MK趋势检验
# 描述：使用Mann-Kendall方法对NDVI时间序列进行趋势检验
# 功能：
#   - 读取栅格格式的NDVI时间序列数据
#   - 对每个像素点进行MK趋势检验
#   - 根据检验结果和Sens斜率对趋势进行分类
#   - 计算全局平均时间序列的累积距平和MK统计量
#   - 生成趋势变化图表
#   - 输出趋势检验结果为栅格文件
# ============================================================================

# ==== 步骤1：导入所需的库 ====
import os  # 文件和路径操作
import numpy as np  # 数组和科学计算
import rasterio  # 栅格数据读写
import pymannkendall as mk  # Mann-Kendall趋势检验
import glob  # 文件模式匹配
import matplotlib.pyplot as plt  # 绘图
from matplotlib import rcParams  # 图表样式设置

# ==== 步骤2：设置中文字体和图表配置 ====
# 设置全局字体为微软雅黑
rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 设置字体为微软雅黑
rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# ==== 步骤3：设置数据路径和获取文件列表 ====
# 设置栅格数据文件夹路径（请根据自己的路径进行修改）
folder_path = r'C:\Users\hp\Pictures\CSDN\MK_Trend\MK_Trend\NDVI'

# 获取所有栅格文件路径
raster_files = glob.glob(os.path.join(folder_path, '*.tif'))

# 按文件名排序，以确保时间顺序正确
raster_files.sort()

# ==== 步骤4：获取栅格数据基本信息 ====
# 读取第一个栅格文件来获取维度信息
with rasterio.open(raster_files[0]) as src:
    rows, cols = src.shape  # 获取行列数
    profile = src.profile  # 获取元数据信息

# ==== 步骤5：创建三维数组并读取时间序列数据 ====
# 创建一个数组来存储所有年份的数据
ndvi_stack = np.zeros((len(raster_files), rows, cols), dtype=np.float32)

# 读取每个栅格文件并存储到数组中
print("开始读取栅格文件...")
for i, raster_file in enumerate(raster_files):
    with rasterio.open(raster_file) as src:
        ndvi_stack[i, :, :] = src.read(1)  # 读取栅格数据
    print(f"已处理 {i + 1}/{len(raster_files)} 个文件 ({((i + 1) / len(raster_files) * 100):.1f}%)")


# ==== 步骤6：定义趋势分类函数 ====
# 定义趋势分类函数，基于MK检验的z值和slope值
def trend_classification_v2(ndvi_stack, rows, cols):
    mk_trend = np.zeros((rows, cols), dtype=np.float32)  # 初始化趋势结果数组
    z_threshold = 1.96  # 对应0.05显著性水平

    total_pixels = rows * cols
    processed_pixels = 0
    print("\n开始进行 Mann-Kendall 趋势检验...")

    # 对每个像素点进行Mann-Kendall检验
    for row in range(rows):
        for col in range(cols):
            pixel_series = ndvi_stack[:, row, col]  # 提取当前像素的时间序列

            # 仅对非NaN值进行检验
            if not np.isnan(pixel_series).all():
                result = mk.original_test(pixel_series)  # 进行MK趋势检验
                z = result.z  # 获取z统计量
                beta = result.slope  # 获取Sens斜率

                # 根据z值和slope值对趋势进行分类
                if beta > 0 and abs(z) > z_threshold:
                    mk_trend[row, col] = 2  # 显著增加
                elif beta > 0 and abs(z) <= z_threshold:
                    mk_trend[row, col] = 1  # 微显著增加
                elif beta == 0:
                    mk_trend[row, col] = 0  # 稳定不变
                elif beta < 0 and abs(z) <= z_threshold:
                    mk_trend[row, col] = -1  # 微显著下降
                elif beta < 0 and abs(z) > z_threshold:
                    mk_trend[row, col] = -2  # 显著下降

            # 进度跟踪
            processed_pixels += 1
            if processed_pixels % (total_pixels // 100) == 0:  # 每完成1%的进度就打印
                progress = (processed_pixels / total_pixels) * 100
                print(f"处理进度: {progress:.1f}%")

    print("Mann-Kendall 趋势检验完成！")
    return mk_trend


# ==== 步骤7：执行趋势分类 ====
# 调用趋势分类函数
mk_trend = trend_classification_v2(ndvi_stack, rows, cols)

# ==== 步骤8：计算全局MK检验统计量 ====
# 计算每个时间点的数据均值（即对所有像元求平均）
data_mean_series = np.nanmean(ndvi_stack, axis=(1, 2))

# 对全局的data均值时间序列进行MK检验
mk_test_result = mk.original_test(data_mean_series)
z_values = mk_test_result.z  # Mann-Kendall 检验的统计量

# ==== 步骤9：计算累积距平 ====
# 计算均值
mean_value = np.nanmean(data_mean_series)

# 计算累积距平值
cumulative_anomaly = np.cumsum(data_mean_series - mean_value)


# ==== 步骤10：计算逐步Mann-Kendall检验统计量 ====
def mk_progressive_test(series):
    """
    计算逐步Mann-Kendall检验统计量

    参数:
        series: 时间序列数据

    返回:
        s: MK统计量序列
    """
    n = len(series)
    s = np.zeros(n)

    for k in range(1, n):
        s[k] = s[k - 1] + np.sign(series[k] - series[:k]).sum()  # 累积计算MK统计量

    return s


mk_statistic_values = mk_progressive_test(data_mean_series)  # 计算MK统计量序列

# ==== 步骤11：绘制图表 ====
plt.figure(figsize=(10, 6))  # 创建图形对象

# 绘制累积距平值
plt.subplot(2, 1, 1)  # 创建子图1
plt.plot(np.arange(1, len(data_mean_series) + 1), cumulative_anomaly, label="累积距平值", color='b')
plt.axhline(0, color='gray', linestyle='--')  # 添加水平线
plt.title('累积距平值')  # 设置标题
plt.xlabel('时间')  # 设置x轴标签
plt.ylabel('累积距平值')  # 设置y轴标签

# 绘制 MK 检验统计量
plt.subplot(2, 1, 2)  # 创建子图2
plt.plot(np.arange(1, len(data_mean_series) + 1), mk_statistic_values, label="MK 检验统计量", color='r')
plt.axhline(0, color='gray', linestyle='--')  # 添加水平线
plt.title('MK 检验统计量')  # 设置标题
plt.xlabel('时间')  # 设置x轴标签
plt.ylabel('MK S值')  # 设置y轴标签

plt.tight_layout()  # 优化布局
plt.show()  # 显示图形

# 保存图像
output_file_img = r'mk_cumulative_anomaly.png'  # 图像输出路径
plt.savefig(output_file_img)  # 保存图像

# ==== 步骤12：保存趋势结果为栅格文件 ====
output_file = r'MK_Trend.tif'  # 栅格输出路径

# 更新数据类型为 float32
profile.update(dtype=rasterio.float32, count=1)  # 设置输出栅格参数

# 写入栅格文件
with rasterio.open(output_file, 'w', **profile) as dst:
    dst.write(mk_trend.astype(rasterio.float32), 1)  # 写入结果数据

print("Mann-Kendall 检验完成，趋势结果已保存至:", output_file)

4 Python代码运行结果