探索Python 融合地学：使用Python一键进行栅格数据Sen+MK长时间序列趋势分析+显著性检验

在长时间序列的植被覆盖NDVI、LAI、气温或降水变化研究中我们经常需要回答两个问题趋势是什么变绿了还是变黄了趋势显著吗是真变了还是偶然波动学术界公认的黄金标准方法便是 Theil-Sen Median (Sen斜率) 结合 Mann-Kendall (MK) 显著性检验。以往我们可能需要借助 MATLAB或者在 ArcGIS 中进行复杂的栅格计算器叠加甚至需要把栅格转成点提取到 Excel 做分析步骤繁琐且容易出错。今天分享一段 Python 代码利用 rasterio 和 numpy 库一键实现从读取多时相 TIFF 影像到输出最终趋势分类图的全过程文末附完整代码和示例数据练习1 原理简介为什么是 Sen MKSens Slope (Sen斜率)对时间序列中所有两两数据点组合逐一计算斜率最终取中位数斜率作为趋势估计值。相比于普通的一元线性回归它对异常值不敏感非常适合处理容易受云雾干扰的遥感数据。其基本公式为判别β 0上升趋势β 0下降趋势β ≈ 0无趋势Mann-Kendall (MK) 检验一种非参数检验方法。它不需要数据服从正态分布能有效检验趋势是否显著即排除随机波动的可能性。统计量S是对时间序列中每对 (xi,xj)ji若后者更大则 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 scipy2完整代码import os import numpy as np import rasterio from scipy.stats import norm from tqdm import tqdm from itertools import combinations # 1. 设置路径 data_dir rC:\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, fLAI_{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), dtypenp.float32) # 3. 读取所有年份数据 print(正在读取数据...) for i, year inenumerate(years): with rasterio.open(os.path.join(data_dir, fLAI_{year}.tif)) as src: data src.read(1).astype(np.float32) # 将原始影像中的 NoData 值统一转换为 np.nan方便后续处理 if src_nodata isnotNone: 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, dtypenp.int8) # 定义插值函数 (处理少量缺失值) deffill_missing_values(y): 简单的线性插值填补 NaN nans np.isnan(y) # 如果没有 NaN直接返回 ifnot np.any(nans): return y # 获取有效值的索引 x np.where(~nans)[0] # 如果有效值太少比如少于 10 年则无法插值返回原数据 iflen(x) 10: return y # 线性插值 y[nans] np.interp(np.where(nans)[0], x, y[~nans]) return y defsens_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 inrange(1, n_years): for k inrange(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 # --- 计算 Sens slope --- sen sens_slope(x) # --- 分类 --- absZ abs(Z) # 阈值判定 if absZ 2.58and sen 0.0005: trend_result[i] 3 elif1.96 absZ 2.58and sen 0.0005: trend_result[i] 2 elif1.645 absZ 1.96and sen 0.0005: trend_result[i] 1 elif absZ 1.645orabs(sen) 0.0005: trend_result[i] 0 # 这里 0 代表不显著是有效结果 elif absZ 2.58and sen -0.0005: trend_result[i] -3 elif1.96 absZ 2.58and sen -0.0005: trend_result[i] -2 elif1.645 absZ 1.96and sen -0.0005: trend_result[i] -1 # 6. 保存 trend_map trend_result.reshape(rows, cols) # 修改点 C: 更新 nodata 为 99 meta.update(dtyperasterio.int8, count1, nodataNODATA_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.012显著增加P 0.051.96 ≤1微显著增加P 0.11.645 ≤0不显著/无变化P ≥ 0.1-1微显著减少P 0.11.645 ≤-2显著减少P 0.051.96 ≤-3极显著减少P 0.01注意代码中设置了一个斜率阈值 0.0005如果 Sen 斜率的绝对值小于该值也会被强制归类为“不显著/无变化”这有助于过滤掉那些微乎其微的噪声变化。好了今天的内容就分享到这里了如果对你有帮助不要忘记了给我们点赞哦