在长时间序列的植被覆盖（NDVI、LAI）、气温或降水变化研究中，我们经常需要回答两个问题：
- 趋势是什么？（变绿了还是变黄了？）
- 趋势显著吗？（是真变了，还是偶然波动？）
学术界公认的黄金标准方法便是 Theil-Sen Median (Sen 斜率) 结合 Mann-Kendall (MK) 显著性检验。
以往我们可能需要借助 MATLAB，或者在 ArcGIS 中进行复杂的栅格计算器叠加，甚至需要把栅格转成点提取到 Excel 做分析，步骤繁琐且容易出错。
今天分享一段 Python 代码，利用 rasterio 和 numpy 库，一键实现从读取多时相 TIFF 影像到输出最终趋势分类图的全过程！
1 原理简介
为什么是 Sen + MK？
- Sen's Slope (Sen 斜率)：对时间序列中所有两两数据点组合，逐一计算斜率，最终取中位数斜率作为趋势估计值。相比于普通的一元线性回归，它对异常值不敏感，非常适合处理容易受云雾干扰的遥感数据。其基本公式为：
判别：β > 0：上升趋势；β < 0：下降趋势；β ≈ 0：无趋势
- Mann-Kendall (MK) 检验：一种非参数检验方法。它不需要数据服从正态分布，能有效检验趋势是否显著（即排除随机波动的可能性）。统计量 S 是对时间序列中每对 (xi,xj)（j>i），若后者更大则 +1，更小则 -1，相等则 0，最后求和。标准化后得 Z 统计量。公式：
判别：在显著性水平 α = 0.05（双尾）下，|Z| > 1.96** 表示趋势显著，|Z| ≤ 1.96 表示趋势不显著。
这里我们借助 Ai，绘制了一张完整的计算流程示意图
2 核心代码实现
本脚本支持读取文件夹内按年份命名的 TIF 影像（如 LAI_2001.tif...你需要修改你自己的文件命名格式），自动逐像元计算，最终输出一张带有分类结果的 GeoTIFF 图像。
（1）环境依赖
pip install numpy rasterio tqdm scipy
（2）完整代码
import os
import numpy as np
import rasterio
from scipy.stats import norm
from tqdm import tqdm
from itertools import combinations
# === 1. 设置路径 ===
data_dir = r"C:\adan\data\ndvi"
output_path = os.path.join(data_dir, "ndvi_trend.tif")
years = list(range(2001, 2021))
n_years = len(years)
# === 2. 读取第一个文件，获取元数据 ===
sample_path = os.path.join(data_dir, f"LAI_{years[0]}.tif")
with rasterio.open(sample_path) as src:
meta = src.meta.copy()
# 读取原始 nodata 值，如果没有则默认为 None
src_nodata = src.nodata
data_stack = np.zeros((n_years, src.height, src.width), dtype=np.float32)
# === 3. 读取所有年份数据 ===
print("正在读取数据...")
for i, year in enumerate(years):
with rasterio.open(os.path.join(data_dir, f"LAI_{year}.tif")) as src:
data = src.read(1).astype(np.float32)
# 将原始影像中的 NoData 值统一转换为 np.nan，方便后续处理
if src_nodata is not None:
data[data == src_nodata] = np.nan
data_stack[i] = data
# === 4. 预处理 ===
rows, cols = data_stack.shape[1:]
pixels = rows * cols
data_series = data_stack.reshape(n_years, pixels).T
# === 修改点 A: 初始化为特殊值（如 99），避免与 0 (不显著) 混淆 ===
NODATA_VAL = 99
trend_result = np.full((pixels,), NODATA_VAL, dtype=np.int8)
# === 定义插值函数 (处理少量缺失值) ===
def fill_missing_values(y):
"""简单的线性插值填补 NaN"""
nans = np.isnan(y)
# 如果没有 NaN，直接返回
if not np.any(nans):
return y
# 获取有效值的索引
x = np.where(~nans)[0]
# 如果有效值太少（比如少于 10 年），则无法插值，返回原数据
if len(x) < 10:
return y
# 线性插值
y[nans] = np.interp(np.where(nans)[0], x, y[~nans])
return y
def sens_slope(x):
slopes = [(x[j] - x[i]) / (j - i) for i, j in combinations(range(len(x)), 2)]
return np.median(slopes)
# === 5. 逐像素分析 ===
for i in tqdm(range(pixels), desc="计算趋势 (含插值)"):
x = data_series[i]
# === 修改点 B: 宽松的过滤条件 ===
# 1. 如果全是 NaN 或 全是 0，跳过
if np.all(np.isnan(x)) or np.all(x == 0):
continue
# 2. 尝试插值填补缺失年份
x = fill_missing_values(x)
# 3. 插值后再次检查，如果仍有 NaN（说明有效年份太少插值失败），则跳过
if np.any(np.isnan(x)):
continue
# --- 计算 Mann-Kendall ---
S = 0
for j in range(1, n_years):
for k in range(j):
if x[j] > x[k]:
S += 1
elif x[j] < x[k]:
S -= 1
VAR_S = n_years * (n_years - 1) * (2 * n_years + 5) / 18
if S > 0:
Z = (S - 1) / np.sqrt(VAR_S)
elif S < 0:
Z = (S + 1) / np.sqrt(VAR_S)
else:
Z = 0
# --- 计算 Sen's slope ---
sen = sens_slope(x)
# --- 分类 ---
absZ = abs(Z)
# 阈值判定
if absZ >= 2.58 and sen > 0.0005:
trend_result[i] = 3
elif 1.96 <= absZ < 2.58 and sen > 0.0005:
trend_result[i] = 2
elif 1.645 <= absZ < 1.96 and sen > 0.0005:
trend_result[i] = 1
elif absZ < 1.645 or abs(sen) < 0.0005:
trend_result[i] = 0
# 这里 0 代表不显著，是有效结果
elif absZ >= 2.58 and sen < -0.0005:
trend_result[i] = -3
elif 1.96 <= absZ < 2.58 and sen < -0.0005:
trend_result[i] = -2
elif 1.645 <= absZ < 1.96 and sen < -0.0005:
trend_result[i] = -1
# === 6. 保存 ===
trend_map = trend_result.reshape(rows, cols)
# === 修改点 C: 更新 nodata 为 99 ===
meta.update(dtype=rasterio.int8, count=1, nodata=NODATA_VAL)
with rasterio.open(output_path, "w", **meta) as dst:
dst.write(trend_map.astype(rasterio.int8), 1)
print("✅ 趋势分析完成！")
3 结果对比
代码运行结束后，你会得到一个 TIF 文件。将其拖入 ArcGIS 或 QGIS 中，根据像元值进行唯一值渲染（Unique Values）。
像元值对应的统计学意义如下表：
| 像元值 (Pixel Value) | 趋势类型 | 显著性水平 | Z 值范围 |
|---|---|---|---|
| 3 | 极显著增加 | P < 0.01 | |
| 2 | 显著增加 | P < 0.05 | 1.96 ≤ |
| 1 | 微显著增加 | P < 0.1 | 1.645 ≤ |
| 0 | 不显著/无变化 | P ≥ 0.1 | |
| -1 | 微显著减少 | P < 0.1 | 1.645 ≤ |
| -2 | 显著减少 | P < 0.05 | 1.96 ≤ |
| -3 | 极显著减少 | P < 0.01 |
注意：代码中设置了一个斜率阈值 0.0005，如果 Sen 斜率的绝对值小于该值，也会被强制归类为'不显著/无变化'，这有助于过滤掉那些微乎其微的噪声变化。
