Sen斜率与Mann-Kendall检验：时间序列稳健趋势分析实战指南

搞地学、遥感、水文、气象数据分析的老哥老姐们，Sen+MK趋势分析这组词，估计十有八九都熟。它基本是“时间序列趋势检测”里的默认组合了：Sen负责估算变化速率，MK负责判断趋势显著性，两个一配合，既能告诉你“涨了还是跌了”，又能告诉你“这个涨跌能不能信”。我自己日常处理的不是NDVI就是降水、土壤湿度，凡是论文里写“分析近二十年演变趋势”的，基本都用这套方法交差。这篇文章就把我从原理理解、数据准备到Python跑批的完整流程，连同踩过的坑一起拆开，给刚开始接触趋势分析的学生，也给每一个被长序列数据折磨的工程师。

先给个结论：这套方法不挑数据分布、对离群值不敏感、实现成本低，属于那种“你不需要懂太多数学也能用对”的工具。但它也有边界，比如对自相关数据、季节数据、短序列数据，直接套用会出大问题。下面我会先把原理讲透，再给完整代码和参数选择思路，最后把高频问题列成排查清单，照着做基本能避掉九成以上的坑。

1. 为什么是“Sen+MK”：这对组合解决什么问题

1.1 常规线性回归的尴尬：被离群值牵着走

拿到一条时间序列，多数人第一反应是做一元线性回归，看斜率是否显著。这在教科书上没问题，真实数据却常常不配合：某一年降水特别大、某一年NDVI因为传感器问题突然跳变、或者源数据本身就带着几个极端噪声点。最小二乘回归对这类离群值非常敏感，一个异常点就足以把斜率从正拽成负，或者从显著拉成不显著。

我印象很深的一次教训：某站点20年降水序列，我做线性回归得到的趋势是“每年增加约3.1毫米，p=0.02”，看上去很漂亮。后来回查原始记录，发现第7年有一次极端特大暴雨，数值比正常年份高出好几倍。把这个点去掉再算，回归斜率直接变成负的。一个点翻转整个结论，这显然不合理。后来又用Sen斜率一算，数字稳定在“每年下降约0.8毫米”，因为Sen拿的是所有点对斜率的中位数，那个极端暴雨点最多影响它参与的少数几个点对，不足以撼动中位数。这也是为什么在真实项目中，我宁可把线性回归当作初筛，最终结论一律以Sen+MK为准。

如果你不理解中位数的意义，可以这样类比：一群人里混进一个“年薪百亿”的，算平均收入瞬间被拉爆，但算收入中位数，大家的生活水平几乎不受影响。Sen斜率就是“收入中位数版”的趋势估计，它不追求完美拟合，只求给一个不被异常值绑架的稳健方向。

1.2 Sen斜率估计：把所有点对斜率排个队

Sen斜率（也叫Theil-Sen估计）计算起来非常朴素：对于一条长度为n的序列，任意取两个时间点i和j（i < j），算出一个斜率：

β = median( (x_j - x_i) / (j - i) ), i < j

就是说，20年的数据共有20×19/2 = 190个点对，每个点对产生一个斜率，最后取这190个斜率的中位数。这样估计出来的β含义很直观：在整段时间内，该指标平均每年变化多少单位。如果时间间隔不是等间隔的，分母要用实际的时间差，这在部分站点数据里会用到。

为什么要强调“中位数”而不是“均值”？因为点对斜率里同样可能混入离群值。取均值的稳健性和最小二乘差不多，取中位数才能把离群点对的影响压到最低。Sen斜率本身不要求数据服从正态分布，不要求等方差，只要求数据之间近似独立，这比线性回归的假设宽松得多。它在线性回归里还有一个变体叫Theil-Sen稳健回归，原理也是拿所有点对斜率的中位数作为整体斜率，本文章节3.2的代码可以直接迁移过去。

1.3 Mann-Kendall检验：只比大小，不看数值

Sen斜率告诉我们变化方向和大小，但没说这个变化“可不可信”。比如20年里降水的Sen斜率是-0.5毫米/年，但这个数字可能纯粹是随机波动凑出来的，统计上不显著。Mann-Kendall检验就是干这件事的。

MK检验属于非参数检验，原假设是“序列不存在单调趋势”。它不关心数值的具体量级，只关心数值之间的相对大小关系。做法是：对序列里所有点对(x_i, x_j)，i < j，比较大小，x_j大于x_i记1分，小于记-1分，相等记0分，所有分数加总得到统计量S：

S = Σ sign(x_j - x_i)

如果S是一个很大的正数，说明后面的值系统性地比前面大，序列整体在上升；如果S负得很厉害，说明在下降。光有S还不够，还要知道S在“完全无趋势”的假设下会怎样波动，于是构造标准化统计量Z：

当S > 0时，Z = (S - 1) / √Var(S) 当S < 0时，Z = (S + 1) / √Var(S) 当S = 0时，Z = 0

其中Var(S) = [n(n-1)(2n+5) - Σ_t t_i(t_i-1)(2t_i+5)] / 18，后一项是“并列值修正”。如果数据里有大量完全相等的点，比如降水经常出现0.0毫米的月份，这个修正必须做，否则方差被低估，p值虚小，容易得出“假显著”结论。

Z值近似服从标准正态分布，所以双侧检验p值就是2倍的尾部概率：p = 2×(1 - Φ(|Z|))。p小于0.05，我们就说趋势显著；小于0.01，说极显著。MK检验不要求数据正态，不怕离群值，对大多数单调趋势的检测都很有效，所以成了环境变化研究里的常青树。

1.4 为什么两者总是成对出现

Sen和MK能经常合体，是因为它们优势互补、方法论基因一致。

Sen算出“每年增加0.003的NDVI”，MK负责告诉你“p=0.42，这个增加很可能是噪声”。反过来，MK说“z=2.31，p=0.02，存在显著趋势”，Sen则告诉你是每年1毫米还是每年50毫米。只用一个都容易误读：一个“显著”但没有量级，一个“有量级”但不可信。在实际项目中，我还会顺带输出MK计算过程中得到的Kendall Tau，它介于-1和1之间，表示趋势的强度。写报告时，表格里放“Sen斜率、p值、Tau值”三列，基本就齐活了。

2. 动手前的准备：数据形态、工具选型与几项关键约定

2.1 数据应该是什么样：一维时间序列和多波段栅格

Sen+MK的输入，按项目场景分两种。第一种是单点时间序列，比如某个气象站的逐年降水量，或者某个采样点的水质指标，一维数组，长度为年份数。第二种是空间栅格，比如MODIS NDVI的年合成产品，20年就是20个波段，每个波段代表一年的空间分布，合成一个多波段GeoTIFF。处理栅格时，本质上是把每个像素对应的“时间维数组”单独拎出来做一次Sen+MK，再把结果写回成一个单波段栅格。

无论哪种形态，有一个前提必须搞清楚：时间分辨率是什么。年度数据用普通MK没问题，月度、周度、日度数据如果直接套普通MK，结果基本会失真。因为普通MK的方差公式假设数据是独立同分布的，而季节数据存在强周期，后面第4.3节会专门讲正确的处理方式。另外，序列里允许有缺失值，但要注意有效样本量，少于一定数量就别勉强出结论。

2.2 Python环境与依赖安装

我现在的标准环境是Python 3.9以上，主要依赖：numpy、scipy、pandas、pymannkendall、rasterio、matplotlib。前三个负责基础计算和数据整理，pymannkendall是趋势检验的核心，rasterio处理栅格读写，matplotlib出图。

pip install numpy scipy pandas pymannkendall rasterio matplotlib

pymannkendall这个库非常完整，原始MK、季节MK、自相关修正MK、预白化MK全都有，函数名也直观。如果你不用Python，R语言里也有zyp和trend包，功能类似，但代码风格不同。这里我只讲Python路径，毕竟栅格批处理这一块Python生态确实更顺手。

2.3 两个常用实现路径：手搓公式与现成函数

做趋势分析最忌讳的就是黑盒使用。我的习惯是：先用现成函数快速出结果，再手写一遍核心公式验证理解。项目紧的时候直接调库，项目松的时候至少把S统计量和Sen斜率自己算一遍，这样出了意外结果才知道往哪个方向排查。

pymannkendall最核心的接口是original_test，它返回一个命名元组，里面包含：trend（"increasing"/"decreasing"/"no trend"）、h（是否拒绝原假设）、p（p值）、z（Z统计量）、Tau（Kendall Tau）、s（S统计量）、var_s（S的方差）、slope（Sen斜率）、intercept（截距）。这里面的slope就是Theil-Sen斜率，截距可以用中位数法另行估计，多数场景下我们不关心截距，只拿斜率。

手写版本虽然慢，但逻辑透明。尤其当你需要处理自定义的“时间间隔不等”或“某些年份缺测”时，现成函数帮不了你，手写公式改起来更灵活。下面第三节我先给手写实现，再给调库版本，两条路都能跑通。

3. 完整实操流程：从一条曲线到上百万像素

3.1 第一步：清洗数据，先给序列做“体检”

任何趋势分析，数据清洗不过关都会白费功夫。拿到序列后我一般按这个顺序检查：

• 时间轴是否升序排列，有没有重复年份。
• 有没有物理上不可能的值，比如降水为负、NDVI超出[-1,1]、温度超过区域历史极值。
• 缺失值如何处理。如果只是零星缺测，可以直接跳过；如果连续缺很多年，比如20年缺了6年，我建议放弃这个站点或把这个像素标记为无数据。
• 是否含有大量完全相同的值。降水序列里一堆0.0毫米是正常的，MK的并列值修正会处理，但如果你用的是别人写的不带tie修正的代码，这里就是个坑。

清洗代码很简单，但作用很大：

import numpy as np import pandas as pd def clean_series(series, min_valid=15): series = np.asarray(series, dtype=float) valid = ~np.isnan(series) if valid.sum() < min_valid: return None # 把NaN直接剔除，注意如果是栅格，这里返回的是有效值序列 return series[valid]

我的经验是：20年数据至少要有15个有效年份才出趋势结果；30年数据至少25个有效年份；少于这个阈值，MK检验的功效会明显下降，结论很容易站不住脚。所谓“检验功效下降”，通俗说就是：趋势可能真实存在，但你手里的数据点太少，无法把它和随机波动区分开。

这里还要提一个反直觉的坑：不要轻易对缺失年份做线性插值。插值会人为拉平数据，降低序列的真实波动，导致Sen斜率被低估，而且插值引入了额外的数据依赖，MK检验的独立性假设会被破坏。只要缺测不是连续多年，宁可跳过缺失点直接算，也别“补”出一个看似完整实则虚假的序列。

3.2 第二步：计算Sen斜率，所有点对挨个算

数据清洗完后，第一步先算Sen斜率。下面是手写版，逻辑和公式一一对应：

import numpy as np def sen_slope(series): series = np.asarray(series, dtype=float) valid = ~np.isnan(series) x = series[valid] n = x.size if n < 3: return np.nan slopes = [] for i in range(n - 1): for j in range(i + 1, n): # 等间隔时间序列，分母就是j-i；非等间隔可改用真实时间差 slopes.append((x[j] - x[i]) / (j - i)) return np.median(slopes)

你可能会担心这个双重循环太慢。确实，n=20时有190个点对，n=50时则上升到1225个，但单条序列毫无压力。真正需要担心的是栅格场景，几百万像素乘几千个点对，那才需要优化，后面第3.4节会讲。

模拟数据验证一下：我造一条20年、真实趋势为每年增加0.08的序列，并在第7年塞进一个4.2的离群值。

rng = np.random.default_rng(42) n = 20 t = np.arange(n) signal = 0.5 + 0.08 * t + rng.normal(0, 0.3, n) signal[7] = 4.2 # 人为离群值 print(f"Sen斜率: {sen_slope(signal):.4f} 单位/年")

输出大概是“Sen斜率: 0.0827 单位/年”，非常接近真实的0.08。如果换成普通最小二乘斜率，很容易被第7年的离群点拉高到0.2以上。这就直观说明了“稳健”两个字的分量。

3.3 第三步：跑Mann-Kendall检验，先手写，再调库

MK检验的手写实现稍微长一点，但核心就是算S、算方差、算Z、算p。我建议你至少完整写一遍，后面不管换什么方法都能知其所以然。

from collections import Counter from scipy import stats import numpy as np def mk_test(series): series = np.asarray(series, dtype=float) valid = ~np.isnan(series) x = series[valid] n = x.size if n < 3: return np.nan, np.nan, np.nan s = 0 for i in range(n - 1): for j in range(i + 1, n): if x[j] > x[i]: s += 1 elif x[j] < x[i]: s -= 1 # 相等贡献0，不做事就是0 # 并列值修正：如果有相同的数，调整方差 counts = Counter(x) tie_term = 0 for v, c in counts.items(): if c > 1: tie_term += c * (c - 1) * (2 * c + 5) var_s = (n * (n - 1) * (2 * n + 5) - tie_term) / 18.0 if s > 0: z = (s - 1) / np.sqrt(var_s) elif s < 0: z = (s + 1) / np.sqrt(var_s) else: z = 0.0 p = 2 * (1 - stats.norm.cdf(abs(z))) return s, z, p s, z, p = mk_test(signal) print(f"S = {s}, Z = {z:.3f}, p = {p:.4f}")

把这个手写结果和pymannkendall的original_test对比一下，两者应该完全一致。这样你就有了两个可信来源，以后不管谁问，都能说清楚统计量是怎么算出来的。

实际项目里，调用库函数更快：

import pymannkendall as mk result = mk.original_test(signal) print(result.trend) # increasing / decreasing / no trend print(result.p) # p值 print(result.z) # Z统计量 print(result.slope) # Theil-Sen斜率 print(result.Tau) # Kendall Tau

original_test返回的slope就是第3.2节手写函数的结果，两者应该对得上。注意，result.trend对“是否拒绝原假设”的判断使用的是0.05显著水平，如果你想要0.1或者0.01，别直接用trend字段，改用p值自己做判断。

3.4 第四步：批量处理栅格，逐像素确实慢，但别怕

真正让Sen+MK项目从“练习”变成“实战”的，是空间栅格批处理。假设你有一个NDVI年合成多波段GeoTIFF，shape是(20, rows, cols)，20个波段代表20年。目标就是为每一个像素计算Sen斜率和MK的p值，输出两张单波段栅格。

最直观的写法是直接用rasterio读成numpy数组，然后双重循环遍历每个像素，把每像素的时间序列取出来跑original_test。代码如下：

import rasterio import numpy as np import pymannkendall as mk with rasterio.open('ndvi_annual_stack.tif') as src: arr = src.read() # shape: (n_year, rows, cols) profile = src.profile profile.update(count=1, dtype='float32', nodata=np.nan) rows, cols = arr.shape[1], arr.shape[2] slope_raster = np.full((rows, cols), np.nan, np.float32) p_raster = np.full((rows, cols), np.nan, np.float32) z_raster = np.full((rows, cols), np.nan, np.float32) min_valid = 15 # 有效样本阈值 for r in range(rows): for c in range(cols): series = arr[:, r, c].astype(float) valid = ~np.isnan(series) if valid.sum() < min_valid: continue res = mk.original_test(series[valid]) slope_raster[r, c] = res.slope p_raster[r, c] = res.p z_raster[r, c] = res.z if (r + 1) % 100 == 0: print(f"processed {r + 1}/{rows} rows")

这段代码能跑，但性能完全取决于影像尺寸。1000×1000的影像，循环里要处理100万个像素，每个像素20年数据约产生190个点对，总共约1.9亿次比较，Python纯循环跑起来挺酸爽。如果影像更大、年份更长，建议做三件事：

一是先做一个有效像元掩膜，只对“至少连续或非连续有效值达到阈值”的像素计算，空值区域直接跳过。二是用rasterio的窗口读取，按块处理，而不是整个读入内存。三是把像素级MK逻辑写成numba加速版本，或者用multiprocessing并行，这在小范围样机上也有效。

我实际处理一个全国范围的30年NDVI栅格（约8000列×6000行）时，纯Python逐像素跑了大半天。后来改成“先生成有效像元mask + 只遍历有效像素 + numba加速核心函数”，时间压缩到十几分钟。所以提醒一句：如果你只跑一次、影像不大，直接用上面的短代码没问题；如果要把这个流程固化到生产环境，一定要上优化手段。

3.5 第五步：结果输出与显著性分级

计算完成后，写回栅格文件：

with rasterio.open('output_slope.tif', 'w', **profile) as dst: dst.write(slope_raster, 1) with rasterio.open('output_pvalue.tif', 'w', **profile) as dst: dst.write(p_raster, 1)

但是，光输出原始slope和p值，最终画图时还要再分类。我常用的分级规则如下：

出图时，我习惯用暖色表示上升、冷色表示下降，不显著区域统一用浅灰或透明。如果你处理的指标有物理范围，比如NDVI理论范围是-1到1，别忘了在色带数值上标注清楚单位，比如“Sen斜率: 0.001/年”，避免读者误读成“每年增加1%”。

另外我还会生成一张“有效样本量”栅格，当作辅助诊断图。当某些区域趋势看上去很奇怪时，第一件事就是查那里的有效年份数是不是太少。很多时候一个显著性“突变区”其实是缺测造成的假象，有了样本量栅格，这种问题一目了然。

4. 问题排查与避坑实录

4.1 有效样本不足和NaN，趋势结果不可信

这是栅格处理中最常见的坑。MODIS NDVI常因为云覆盖导致某些年份某些像素为无效值；降水数据遇到传感器故障也可能缺测。如果缺失年份太多，MK检验的方差估计会失真。举例来说，20年数据里只有8年有效，序列长度和“自由度”都不足，即使真的存在趋势，也很难检测出来。

我的做法是给每个像素设一个阈值，比如有效样本≥15才返回结果，低于阈值直接输出NaN并mask掉。这个阈值不是拍脑袋，而是从检验功效角度考虑的：n=15时MK对中等强度的趋势尚能检测，n=8时基本什么都检不出来。如果你的时间序列本身就是短序列，比如只有10年，那至少也要有8个有效点才勉强能看，而且结论要在文中说明“受限于样本量，结论仅供参考”。

还有一点：即使有效样本数达标，也不能把缺测年份的数值填成0或者均值，这会人为扭曲序列分布。宁可用有效值序列，也不要补出一个伪完整序列。

4.2 时间自相关导致MK误判，假阳性问题

水文气候序列里经常存在“年际持续性”，比如连续几年偏暖、连续几年偏干。这种自相关会让MK检验的方差被低估，Z值虚高，p值虚小，结果就是“明明没有趋势，却检出了显著趋势”，俗称假阳性。这个问题在有滞后一年自相关的序列里尤其明显。

怎么判断？先算一阶自相关：

def lag1_corr(x): valid = ~np.isnan(x) y = x[valid] if y.size < 5: return np.nan return np.corrcoef(y[:-1], y[1:])[0, 1]

如果算出来lag1大于0.2到0.3，我通常会改用pymannkendall里带自相关修正的方法。这个库提供了几个现成接口：

我的经验阈值是：lag1 > 0.2就使用trend_free_prewhitening_modification_test，否则用original_test。这个判断不算严谨，严谨做法是对lag1做显著性检验，但实际业务的批量处理中，我优先选“不轻易漏掉趋势，也不轻易误报趋势”的平衡方案。需要说明的是，预白化也有副作用，它可能削弱真实的低频趋势信号，所以如果原始数据本身有很强的趋势，先预白化再检验可能会让p值偏大。这也是为什么“任何统计方法都不能无脑套”的原因。

4.3 季节数据不能用普通MK，直接用Seasonal MK

月度、周度、日度数据如果带明显季节性，用普通MK会出大问题。最直观的例子：某地温度序列1月-10度、7月30度，这种季节波动幅度远大于年际趋势，MK比较点对大小关系时，季节差异会淹没趋势信号，统计量方差被季节周期污染，结果可能既不稳定也不可信。

正确做法是使用季节性MK。它的核心是：先按“季节/月份”把序列拆开，同一季节内正常计算MK的S统计量，然后把所有季节的S统计量合并，综合计算趋势。pymannkendall提供现成接口：

# monthly_series是已经按时间顺序排好的逐月数据，n_season表示周期长度，月数据就是12 result = mk.seasonal_test(monthly_series, 12)

同理，季度数据用period=4，周数据用period=52。Seasonal MK的优势是既保留了季节内的信息，又不会被年内的季节波动干扰。如果你手头的月度数据只是“用来分析长期趋势”，另一个更朴素但实用的方案是：把每年同一个月的值单独拿出来，组成12条年度序列分别做普通MK，再从整体上综合判断。这个方法管道简单，但结果不如Seasonal MK统一。

还有一个隐蔽坑：如果你的数据是“生长季NDVI”这类每月或每旬合成的，那么不同年之间“第几期”的对应关系必须严格一致，否则Seasonal MK的分组逻辑就会错位。

4.4 栅格计算慢、内存爆，几个优化技巧

跑大范围栅格时，我踩过的坑大体分两类：内存不足和计算过慢。内存问题多半来自直接把整个多波段tif读进内存。一个10000×10000×30的float32影像，理论内存占用约12GB，普通电脑直接爆掉。解决办法很简单：用rasterio的窗口读取，分块循环：

import rasterio from rasterio.windows import Window from rasterio.transform import Affine block_size = 512 with rasterio.open('ndvi_annual_stack.tif') as src: rows, cols = src.height, src.width profile = src.profile profile.update(count=1, dtype='float32', nodata=np.nan) for row0 in range(0, rows, block_size): for col0 in range(0, cols, block_size): window = Window(col0, row0, min(block_size, cols - col0), min(block_size, rows - row0)) block = src.read(window=window) # shape: (n_time, block_h, block_w) # 对block做像素级Sen+MK，再把结果写回对应窗口

计算过慢则取决于你点对循环的效率。如果年份数超过40，每个像素点对数量至少780对，全国范围内就是几十亿次比较。我的建议是：核心函数用numba的njit装饰，或者用numpy的向量化矩阵构建点对差分（以牺牲部分内存换速度），再或者用multiprocessing多进程并行。对单个项目而言，先用几万像素小范围测试正确性，再上全图，避免程序跑了两小时才发现结果全是NaN。

4.5 p值显示和解读，别掉进“科学计数法”的坑

MK输出的p值有时候是1e-16这种极小数字，很多人偷懒写成“p = 0.000”。这个习惯在正式报告里其实很危险，因为“p=0”意味着“绝对不可能发生”，这在统计上不成立。正确写法是“p < 0.001”或“p = 2.3e-16”。

更重要的是理解“统计显著”不等于“实际显著”。一个20年NDVI趋势，Sen斜率只有每年0.0005，p值却高达0.001，统计上显著，因为序列非常平滑、噪声极小。但这个量级在生态学上可能完全没有意义。我每次写结论都会强制要求自己把“斜率数值的单位”和“实际物理意义”同时写出来。比如写成“NDVI平均每年上升0.002，约等于20年累计0.04，占研究区多年平均NDVI的2.1%”，这样读者才有一个直观尺度。否则一个苍白的小数点数字，既无法支撑结论，也无法横向对比。

5. 我的实操体会与几条进阶建议

这套流程我用了很多年，最大的感受是：Sen+MK不是“一键出图”的黑盒，它的门槛在原理理解、数据清洗和情境判断。同一份数据，年度序列和月度序列用错方法，结论可能完全相反；同一份序列，是否做自相关修正，p值可能从0.04跳到0.20。所以每次跑完结果，我都建议你花十分钟做三件事：看一眼Sen斜率是否在物理合理范围内，看一眼没通过显著性的区域是不是样本量特别少，再看一眼空间分布是否出现了“地图上突兀的点状异常”——如果某个像素和周围大片区域趋势明显不一致，先别急着下结论，回到原始数据去查那个点是不是掩膜没做干净。

如果数据里明显存在突变，比如某一年之后趋势方向改变了，我会在Sen+MK之外再叠加一个Pettitt突变检验或者分段Sen+MK，把“整体趋势”拆成“前半段”和“后半段”分别描述。这能很好地弥补MK只检测单调趋势的短板，尤其在分析政策干预、极端气候事件影响时特别有用。还有一点，如果你处理的是长时间序列且要做空间对比，建议统一设置随机种子以保证可复现；如果再配上有效的样本量栅格和p值栅格，你的论文审稿人或项目验收方基本找不到“统计方法交代不清楚”的指责点。

最后直接给你一条经验：趋势分析只是第一步，结论要落在“为什么变”上。Sen+MK告诉你哪里变了、变了多少、变化可不可信，但它不告诉你原因。想做归因，还需要结合土地利用、气候因子、人类活动等空间数据做进一步分析，那才是另一个有意思的开始。