Python 栅格数据 Theil-Sen 趋势分析与 Mann-Kendall 显著性检验

遥感时间序列趋势分析：Theil-Sen 与 Mann-Kendall

在气候变化、环境监测及生态演变研究中，准确识别长时间序列数据的趋势至关重要。Theil-Sen Median 趋势分析（Sen 分析）结合 Mann-Kendall 显著性检验（MK 检验）是处理此类问题的经典非参数方法组合。相比传统最小二乘法，该方法对异常值和缺失值具有更强的鲁棒性，且无需数据服从正态分布。

原理概述

1. Theil-Sen Median 趋势分析

Sen 分析通过计算所有点对斜率的中位数来评估趋势。对于时间序列 $ET_1, ET_2, ..., ET_n$，任意两点 $(i, j)$ 的斜率计算公式为：

其中 $ET_i$ 和 $ET_j$ 分别为不同时间点的数据值，$Q$ 为中位数斜率。若 $Q > 0$ 表示上升趋势，$Q < 0$ 表示下降趋势。

2. Mann-Kendall 显著性检验

MK 检验用于判断 Sen 分析得出的趋势是否具有统计学意义。它通过比较时间序列中每对数据点的符号来计算统计量 $S$：

其中 $sgn(ET_j - ET_i)$ 为符号函数。根据 $S$ 及其方差可进一步计算 $Z$ 值：

3. 结果判读

在 0.05 置信水平下，依据 $Z$ 值判断显著性：

• 显著上升：$Z > 1.96$
• 显著下降：$Z < -1.96$
• 不显著：$-1.96 \le Z \le 1.96$

代码实现

以下 Python 脚本实现了从数据加载、Sen 斜率计算、MK 检验到结果重分类的全流程。代码支持逐像元进度可视化，并针对空值进行了处理。

准备工作

确保已安装 rasterio, numpy, tqdm, matplotlib 等依赖库。数据文件夹结构建议如下：

核心代码

import os
import rasterio
import numpy as np
from tqdm import tqdm
import matplotlib.pyplot as plt
from pathlib import Path

# 获取当前工作目录
base_path = os.getcwd()

# 定义数据和结果路径
data_path = os.path.join(base_path, 'data')
result_path = os.path.join(base_path, 'results')

# 创建结果目录（如果不存在）
os.makedirs(result_path, exist_ok=True)

# 定义年份范围
start_year = 2000
end_year = 2020
cd = end_year - start_year + 1  # 时间跨度

try:
    # 读取第一个栅格文件以获取元数据和尺寸信息
    first_file = os.path.join(data_path, f'{start_year}.tif')
    with rasterio.open(first_file) as src:
        a = src.read(1)  # 读取栅格数据
        transform = src.transform  # 栅格的空间转换信息
        metadata = src.meta.copy()  # 获取元数据
        
        # 修复：更新元数据中的数据类型为 float32，解决 dtype 不匹配问题
        metadata.update({
            'dtype': 'float32',
            'nodata': np.nan  # 设置空值为 NaN
        })
        
        # 获取栅格尺寸
        m, n = a.shape
        print(f"栅格尺寸：{m} x {n}")

    # 创建空数组来存储所有年份的数据
    # 使用三维数组以避免展平和重塑操作
    all_data = np.full((m, n, cd), np.nan)

    # 加载每一年的数据
    print("正在加载每一年的栅格数据...")
    for i, year in enumerate(tqdm(range(start_year, end_year + 1), desc="加载年度数据")):
        filename = os.path.join(data_path, f'{year}.tif')
        with rasterio.open(filename) as src:
            all_data[:, :, i] = src.read(1)

    # 创建结果数组
    sen_result = np.full((m, n), np.nan)
    valid_pixels = 0
    
    # 优化：仅处理有效像素（至少有一个非空值的像素）
    valid_mask = np.any(all_data > 0, axis=2)
    total_valid = np.sum(valid_mask)

    # 计算 Sen's Slope 趋势
    # 这是计算最耗时的部分，使用 tqdm 显示进度
    print("正在计算 Sen's Slope 趋势...")
    for i in tqdm(range(m), desc="计算 Sen's Slope"):
        for j in range(n):
            if valid_mask[i, j]:
                data = all_data[i, j, :]
                # 检查数据是否包含有效值
                if np.all(~np.isnan(data)) and np.min(data) > 0:
                    slopes = []
                    for k1 in range(1, cd):
                        for k2 in range(k1):
                            # 计算变化率
                            slope = (data[k1] - data[k2]) / (k1 - k2)
                            slopes.append(slope)
                    # 取中位数作为最终斜率值
                    if slopes:
                        sen_result[i, j] = np.median(slopes)
                        valid_pixels += 1

    print(f"有效像素数：{valid_pixels}/{total_valid} (共{m*n}像素)")

    # 设置输出路径并将 Sen's Slope 结果保存为 GeoTIFF
    sen_output_path = os.path.join(result_path, '基于 sen 的 ET 变化趋势.tif')
    
    # 确保结果为 float32 类型
    sen_result = sen_result.astype('float32')
    
    # 保存 Sen's Slope 结果
    with rasterio.open(sen_output_path, 'w', **metadata) as dst:
        dst.write(sen_result, 1)
    print('Sen\'s Slope 处理完成！结果已保存至:', sen_output_path)

    # 计算 Mann-Kendall 检验
    print("正在计算 Mann-Kendall 检验...")
    mk_result = np.full((m, n), np.nan)
    
    for i in tqdm(range(m), desc="计算 Mann-Kendall 检验"):
        for j in range(n):
            if valid_mask[i, j]:
                data = all_data[i, j, :]
                if np.all(~np.isnan(data)) and np.min(data) > 0:
                    sgnsum = []
                    for k1 in range(1, cd):
                        for k2 in range(k1):
                            # 计算符号差异
                            sgn = np.sign(data[k1] - data[k2])
                            sgnsum.append(sgn)
                    # 计算符号差异的总和
                    mk_result[i, j] = np.sum(sgnsum)

    # 计算 Z 值
    print("正在计算 Z 值...")
    vars_mk = cd * (cd - 1) * (2 * cd + 5) / 18
    z_scores = np.full((m, n), np.nan)
    
    # 处理不同情况的 Z 值计算
    z_scores[~np.isnan(mk_result) & (mk_result == 0)] = 0
    z_scores[~np.isnan(mk_result) & (mk_result > 0)] = (mk_result[~np.isnan(mk_result) & (mk_result > 0)] - 1) / np.sqrt(vars_mk)
    z_scores[~np.isnan(mk_result) & (mk_result < 0)] = (mk_result[~np.isnan(mk_result) & (mk_result < 0)] + 1) / np.sqrt(vars_mk)

    # 保存 Mann-Kendall 检验结果
    mk_output_path = os.path.join(result_path, 'MK 检验结果.tif')
    z_scores = z_scores.astype('float32')
    
    with rasterio.open(mk_output_path, 'w', **metadata) as dst:
        dst.write(z_scores, 1)
    print('Mann-Kendall 检验处理完成！结果已保存至:', mk_output_path)

    # 对 Sen's slope 和 MK 检验结果进行重分类
    print("正在进行重分类...")
    S2 = np.full((m, n), np.nan)
    M2 = np.full((m, n), np.nan)
    
    # 重分类 Sen's slope 结果
    S2[np.isnan(sen_result)] = -9999
    S2[~np.isnan(sen_result) & (sen_result <= -0.0005)] = -1
    S2[~np.isnan(sen_result) & (sen_result >= 0.0005)] = 1
    S2[~np.isnan(sen_result) & (sen_result > -0.0005) & (sen_result < 0.0005)] = 0
    
    # 重分类 MK 检验结果
    M2[z_scores > 1.96] = 2
    M2[~np.isnan(z_scores) & (z_scores <= 1.96)] = 1
    
    # 计算最终重分类结果
    reclassify = (S2 * M2).astype(np.int16)
    
    # 设置输出路径并将重分类结果保存为 GeoTIFF
    reclass_output_path = os.path.join(result_path, '重分类结果.tif')
    
    # 更新元数据中的数据类型为 int16
    metadata.update({
        'dtype': 'int16',
        'nodata': -9999
    })
    
    # 保存重分类结果
    with rasterio.open(reclass_output_path, 'w', **metadata) as dst:
        dst.write(reclassify, 1)
    print('重分类处理完成！结果已保存至:', reclass_output_path)

except Exception as e:
    print(f"处理过程中发生错误：{str(e)}")

运行效果

程序执行过程中会显示详细的进度条，方便监控长耗时任务。最终生成的栅格文件可直接在 GIS 软件中查看。

数据说明

本文示例使用了 MOD16A2H 2000-2020 年宁夏部分区域的蒸散发数据。实际使用时，请确保输入数据的时间序列完整且格式一致（如 GeoTIFF）。