Python | 栅格数据Sen+MK长时间序列趋势分析+显著性检验保姆级教程

Ne0inhk

23 Mar 2026 — 10 min read

遥感代码星球的第002篇代码分享

今天给大家分享的是:

Theil-Sen Median趋势分析

Mann-Kendall显著性检验

基于遥感数据的植被碳水利用效率时空变化和归因分析[J].生态学报,2024,44(01):377-391.DOI:10.20103/j.stxb.202207232112.

在研究气候变化、环境变化、生态变化等领域时，如何准确分析时间序列数据的趋势是一个重要的课题。Theil-Sen Median趋势分析（简称Sen分析）和Mann-Kendall显著性检验（简称MK检验）是两种常用的非参数方法，它们在处理气候、环境、生态等领域的趋势分析中有着广泛的应用。本文将详细介绍这两种方法的基本原理，并探讨它们如何结合使用以提高趋势分析的准确性。同时结尾附完整代码及案例数据供大家学习使用。

01 原理及公式

1. Theil-Sen Median 趋势分析（Sen分析）

Theil-Sen Median 趋势分析是一种稳健的非参数统计方法，适用于具有异常值或缺失值的时间序列数据。该方法通过计算数据的中位数斜率（即Sen斜率）来评估时间序列的趋势。不同于传统的最小二乘法，Sen分析通过中位数的方式计算斜率，使其对异常值和极端值不敏感，从而能够更稳健地捕捉数据中的长期趋势。

具体的计算方法为(本文以ET为例)：

其中，ETi和ETj分别为时间序列中不同时间点的数据值Sen为中位数斜率。

若Sen>0，则表示数据存在上升趋势；若Sen<0，则表示数据存在下降趋势。

2. Mann-Kendall 显著性检验（MK检验）

Mann-Kendall检验是一种用于检测时间序列趋势的非参数统计方法。与Sen分析不同，MK检验的目的是判断时间序列的趋势是否显著，它不要求数据服从任何特定的分布。MK检验通过对时间序列中每对数据点进行比较，计算出趋势的符号，并利用统计量评估数据的趋势方向和显著性。

MK检验的基本步骤如下：

对时间序列中所有的每对数据点进行比较，得到每对数据的符号（升序为1，降序为-1，相等为0）；

MK检验的检验统计量S反映了时间序列的总体趋势。根据S值的大小及其方差，可以进一步计算出Z值，用于判断ET变化的显著性。计算S的方差，进一步得到Z值，用于判断趋势的显著性：

Z值越大，表示趋势越显著。

计算符号函数的累加和，得到统计量S：

其中，sgn(ETj−ETi)为符号函数。

3. MK显著性检验结果的划分

在进行MK检验时，Z值的大小决定了趋势的显著性水平。在常见的0.05置信水平下，MK检验的显著性结果可以通过Z值来判断：

显著变化: 上升趋势（Z>1.96）或下降趋势（Z<−1.96）
不显著变化：（−1.96≤Z≤1.96）

4. Sen+MK方法的优势与应用

稳健性强：Sen分析能够有效处理含有异常值或缺失值的时间序列数据，适用于多种复杂环境下的数据分析。
显著性检验：MK检验能够评估趋势的显著性，确保得到的趋势结果是有统计学意义的。
适用范围广泛：Sen+MK方法在气候变化、环境监测、生态变化等领域中广泛应用，尤其适用于无法满足正态分布假设的时间序列数据。

02 代码复现

1. 代码可实现：

✅逐像元栅格进度可视化

✅只需要改动标注有✏️的位置

✅可1比1复刻

2. 文件夹结构：

3.具体代码：

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() # 获取元数据，包括GeoTIFFTags # 修复：更新元数据中的数据类型为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) # 使用Sen's Slope方法计算趋势 # 这是计算最耗时的部分，使用tqdm显示进度 print("正在计算Sen's Slope趋势...") valid_pixels = 0 # 优化：仅处理有效像素（至少有一个非空值的像素） valid_mask = np.any(all_data > 0, axis=2) total_valid = np.sum(valid_mask) 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: # 确保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)}")