在气候变化、生态环境等研究中，如何准确分析长时间序列数据的趋势至关重要。Theil-Sen Median 趋势分析（Sen 分析）和 Mann-Kendall 显著性检验（MK 检验）是两种常用的非参数方法，特别适合处理含有异常值或缺失值的遥感数据。本文将梳理这两种方法的原理，并提供完整的 Python 代码实现。
核心方法解析
1. Theil-Sen Median 趋势分析（Sen 分析）
这是一种稳健的非参数统计方法。不同于传统的最小二乘法，Sen 分析通过计算数据的中位数斜率来评估时间序列的趋势，对异常值和极端值不敏感。
具体计算以蒸散发（ET）为例：
其中，ETi 和 ETj 分别为时间序列中不同时间点的数据值，Sen 为中位数斜率。若 Sen > 0，表示上升趋势；若 Sen < 0，表示下降趋势。
2. Mann-Kendall 显著性检验（MK 检验）
MK 检验用于判断时间序列的趋势是否显著，不要求数据服从特定分布。它通过对每对数据点进行比较，计算趋势的符号，并利用统计量评估显著性。
基本步骤如下：
- 比较所有数据点对，得到符号（升序为 1，降序为 -1，相等为 0）；
- 计算统计量 S，反映总体趋势；
- 根据 S 及其方差计算 Z 值，判断显著性。
Z 值越大，趋势越显著。符号函数累加和得到统计量 S：
3. MK 显著性结果划分
在 0.05 置信水平下，依据 Z 值判断：
- 显著变化：上升趋势（Z > 1.96）或下降趋势（Z < -1.96）
- 不显著变化：-1.96 ≤ Z ≤ 1.96
4. 方法优势
- 稳健性强：有效处理含异常值或缺失值的数据。
- 显著性检验：确保趋势结果具有统计学意义。
- 适用范围广：适用于无法满足正态分布假设的时间序列数据。
代码实现
以下代码实现了逐像元栅格进度可视化，只需修改标注位置即可运行。代码逻辑上优化了空值处理和内存占用，适合长时间序列分析。
import os
rasterio
numpy np
tqdm tqdm
matplotlib.pyplot plt
pathlib Path
base_path = os.getcwd()
data_path = os.path.join(base_path, )
result_path = os.path.join(base_path, )
os.makedirs(result_path, exist_ok=)
start_year =
end_year =
cd = end_year - start_year +
:
first_file = os.path.join(data_path, )
rasterio.(first_file) src:
a = src.read()
transform = src.transform
metadata = src.meta.copy()
metadata.update({
: ,
: np.nan
})
m, n = a.shape
()
all_data = np.full((m, n, cd), np.nan)
()
i, year (tqdm((start_year, end_year + ), desc=)):
filename = os.path.join(data_path, )
rasterio.(filename) src:
all_data[:, :, i] = src.read()
sen_result = np.full((m, n), np.nan)
()
valid_pixels =
valid_mask = np.(all_data > , axis=)
total_valid = np.(valid_mask)
i tqdm((m), desc=):
j (n):
valid_mask[i, j]:
data = all_data[i, j, :]
np.(~np.isnan(data)) np.(data) > :
slopes = []
k1 (, cd):
k2 (k1):
slope = (data[k1] - data[k2]) / (k1 - k2)
slopes.append(slope)
slopes:
sen_result[i, j] = np.median(slopes)
valid_pixels +=
()
sen_output_path = os.path.join(result_path, )
sen_result = sen_result.astype()
rasterio.(sen_output_path, , **metadata) dst:
dst.write(sen_result, )
(, sen_output_path)
()
mk_result = np.full((m, n), np.nan)
i tqdm((m), desc=):
j (n):
valid_mask[i, j]:
data = all_data[i, j, :]
np.(~np.isnan(data)) np.(data) > :
sgnsum = []
k1 (, cd):
k2 (k1):
sgn = np.sign(data[k1] - data[k2])
sgnsum.append(sgn)
mk_result[i, j] = np.(sgnsum)
()
vars_mk = cd * (cd - ) * ( * cd + ) /
z_scores = np.full((m, n), np.nan)
z_scores[~np.isnan(mk_result) & (mk_result == )] =
z_scores[~np.isnan(mk_result) & (mk_result > )] = (mk_result[~np.isnan(mk_result) & (mk_result > )] - ) / np.sqrt(vars_mk)
z_scores[~np.isnan(mk_result) & (mk_result < )] = (mk_result[~np.isnan(mk_result) & (mk_result < )] + ) / np.sqrt(vars_mk)
mk_output_path = os.path.join(result_path, )
z_scores = z_scores.astype()
rasterio.(mk_output_path, , **metadata) dst:
dst.write(z_scores, )
(, mk_output_path)
()
S2 = np.full((m, n), np.nan)
M2 = np.full((m, n), np.nan)
S2[np.isnan(sen_result)] = -
S2[~np.isnan(sen_result) & (sen_result <= -)] = -
S2[~np.isnan(sen_result) & (sen_result >= )] =
S2[~np.isnan(sen_result) & (sen_result > -) & (sen_result < )] =
M2[z_scores > ] =
M2[~np.isnan(z_scores) & (z_scores <= )] =
reclassify = (S2 * M2).astype(np.int16)
reclass_output_path = os.path.join(result_path, )
metadata.update({
: ,
: -
})
rasterio.(reclass_output_path, , **metadata) dst:
dst.write(reclassify, )
(, reclass_output_path)
Exception e:
()
