每次拿到降水、径流、气温这类长序列数据我第一个想跑的分析大概率都是趋势检验。而在水文、气象、环境领域干了这些年我发现折腾到最后几乎所有人都会回到同一个组合——Sen斜率估计加上Mann-Kendall检验也就是大家常说的SenMK趋势分析。这套组合在我手头处理过几百组站点数据算得上一对久经考验的老搭档今天就把它们彻底拆开讲清楚。这套方法解决的痛点很具体我们拿到一份30年降水序列想知道它到底有没有逐年增加或减少的倾向变化速率又是多少。线性回归也能给答案但回归对异常值太敏感数据又不一定满足正态分布算出来的p值经常不可靠。SenMK恰好绕开了这些限制不要求数据正态、不容易被极端值带偏而且从1960年代提出到现在在各种学科的论文里被反复验证说服力很强。如果你是做气候变化、水文分析、环境质量评估或者农业气象相关工作的这篇内容可以直接帮你把整个分析流程跑通。我会从原理讲到手写代码再讲实操里的坑最后给出一整套可以直接抄作业的案例。1. 为什么趋势分析绕不开SenMK这对组合1.1 和线性回归相比SenMK到底强在哪先说线性回归做趋势分析的缺陷。最小二乘法那条直线本质上是在最小化所有点到直线距离的平方和。这个概念放在正态分布、方差齐性的数据上没问题可水文气象数据偏偏经常带几个异常年份——比如某年特大洪水、极端干旱。一个离群点就会把回归线的斜率拉偏更麻烦的是异常点容易制造伪显著的p值。我遇到过最典型的一个案例某站降雨数据1998年出现一个极端峰值线性回归做出来slope -3.2 mm/年p 0.03看着趋势显著下降。可把这个峰值剔除再算斜率立刻变成-0.8 mm/年p值变成0.41。一条极端年份的数据就能让结论反转这显然不够稳健。Sen斜率估计的处理方式完全不同。它不拟合任何直线而是把所有数据两两组合计算每一对点之间的斜率最后取这些斜率的中位数。中位数的特性决定了它天然抗异常值就算数据里有几个极端年份只要不超过一半对中位数的影响就极其有限。这就是为什么在长时间序列趋势分析里Sen斜率比最小二乘斜率更被认可。MK检验则是跟Sen斜率配套的显著性检验。它不像t检验那样要求数据来自正态总体而是纯粹基于数据的秩排序位置来判断序列是否有单调上升或下降的趋势。两个方法一个估计幅度、一个检验显著性搭配使用刚好补全了对方的短板。1.2 哪些场景最适合这套组合以我的实际经验来看下面这几类数据用SenMK效果最好水文序列径流量、降雨量、蒸发量、水位、泥沙含量气象序列气温、降水、风速、日照时数、极端气候指数环境监测河流水质指标COD、氨氮、总磷、PM2.5浓度、地下水埋深遥感反演产品NDVI、植被覆盖度、土壤湿度、雪水当量这些数据的共同特点是观测时间跨度长一般要求至少10年以上建议20年以上、存在周期性波动、偶尔有缺失值或异常值、分布往往偏态。这些条件凑在一起参数方法基本都会失效非参数方法就成了唯一靠谱的选择。需要注意SenMK分析的是单调趋势不是周期变化。如果你的数据有明显的季节节律比如月径流数据冬天低夏天高那必须先做季节分解或者按相同月份分别分析否则结果会被季节波动污染。后面我会专门讲这个坑。2. Sen斜率与MK检验的原理用手算也能跑通2.1 Sen斜率的计算逻辑和背后的直觉Sen斜率的正式名称是Theil-Sen估计量Theil 1950年提出Sen 1968年推广。假设我们有n个时间点的观测值(x₁, y₁), (x₂, y₂), …, (xₙ, yₙ)其中x通常是时间序号第1年、第2年……y是观测值。计算分两步第一步把所有i j的数据对配对计算斜率slope_ij (y_j - y_i) / (x_j - x_i)这样总共会产生n(n-1)/2个斜率值。比如30年数据就是30×29/2 435个斜率。第二步对这435个斜率取中位数Sens slope median(slope_ij)这个中位数斜率的意义是整个时间段内变量每单位时间通常是每年的典型变化量。比如计算结果是0.35 mm/年意思就是降水量以每年0.35毫米的速率增加。为什么中位数比均值靠谱举个例子一组数据[1, 2, 3, 4, 100]均值是22中位数是3。100这个异常值把均值拉得面目全非但中位数纹丝不动。Sen斜率继承了同样的抗干扰能力单对极端值根本无法主宰最终结果。实际应用里还有一个加分项Sen斜率天然给出了单位方便直接写入报告。不像相关系数只是一个无量纲的数字Slope 3.2 mm/年可以直接说该站降水量以每年3.2毫米的速率减少非常直观。2.2 MK检验的统计量构造过程Mann-Kendall检验Mann 1945, Kendall 1975的核心是构造一个S统计量。对于序列x₁, x₂, …, xₙS的计算公式是S Σ(ij) sign(x_j - x_i)其中sign是符号函数sign(t) 1, t 0 sign(t) 0, t 0 sign(t) -1, t 0换句话说把所有后面的观测值和前面的比一遍后面比前面大加1分后面比前面小减1分相等不得分。把所有分数加总就是S。如果S是一个很大的正数说明后面数据整体比前面大——序列在上升如果S是很大的负数说明序列在下降如果S接近0说明数据随时间没有明确方向波动为主。接下来的关键问题是S大到哪里才算显著这需要用统计推断。在原假设无趋势下S的期望是0而方差Var(S)有明确的解析表达式Var(S) [n(n-1)(2n5) - Σ(t_p)(t_p-1)(2t_p5)] / 18公式里的t_p是数据中第p个并列组的观测个数。如果数据里没有重复值后面的求和项就是0公式退化为Var(S) n(n-1)(2n5) / 18然后用S和Var(S)构造标准化统计量ZZ (S - 1) / sqrt(Var(S)), S 0 Z 0, S 0 Z (S 1) / sqrt(Var(S)), S 0分子里的±1是连续性修正continuity correction目的是让离散的S分布更好地逼近连续的正态分布。这个修正量虽然小但在样本量n不大时会影响最终p值结果很多只抄代码不看原理的人会忽略这一项。Z统计量在大样本下近似服从标准正态分布N(0,1)。判断显著性时只要看Z是否超出标准正态分布的临界值在95%置信水平下|Z| 1.96就拒绝原假设认为趋势显著在90%置信水平下临界值变成1.645。2.3 显著性判读不要只会看p值很多论文汇报MK检验结果时只写p 0.05趋势显著这其实浪费了很多信息。一套完整的判读逻辑应该是方向Z 0表示上升趋势Z 0表示下降趋势。这个方向应该和Sen斜率符号一致如果不一致说明数据有问题需要检查。强度|Z|越大趋势越强。但注意统计上显著不等于实际变化幅度大。30年数据如果每年只增加0.01 mm虽然30年累计0.3 mm在统计上可能显著但实际影响微乎其微。这时候Sen斜率的数值才是真正要重点看的结果。趋势突变MK检验假设趋势是单调的。如果数据先上升后下降比如气温在1990年之前降温、之后升温MK检验的S正负会互相抵消很容易得出无显著趋势的错误结论。遇到这种数据建议先看一眼时间序列图或者累积距平曲线别直接套检验。3. 用Python手写SenMK计算全过程市面上的现成工具包不少但我还是建议自己写一遍核心逻辑这样既能加深理解也方便在特殊情况下改造。下面这个实现不依赖任何水文专用包只要numpy、pandas、scipy就可以跑通。3.1 环境准备与数据组织我用的是Python 3.9以上版本需要的库pip install numpy pandas scipy matplotlib数据格式长这样一个时间列一个数值列year, value 1990, 512.3 1991, 489.7 1992, 531.2 ... 2022, 468.9读进来之后第一件事是排序。按时间升序排列然后检查是否有重复年份。重复年份在MK检验里会造成并列组ties虽然公式里已经考虑了修正但最好还是从源头清理干净。import numpy as np import pandas as pd from scipy import stats def load_data(filepath): df pd.read_csv(filepath) df df.sort_values(year).reset_index(dropTrue) # 检查缺失值和重复年份 assert df[year].is_unique, 存在重复年份请检查数据 assert df[value].notna().all(), 存在缺失值请先处理 return df[year].values, df[value].values year, value load_data(station_precip.csv)3.2 核心函数Sen斜率计算Sen斜率的计算用嵌套循环或者numpy广播都可以。数据量几千个点以下直接双重循环速度也能接受。我为了可读性先用通俗写法后面再给一个向量化版本。def sen_slope(x, y): 计算Theil-Sen斜率 x: 时间序列通常是年份 y: 观测值序列 返回: Sen斜率估计值 n len(x) slopes [] for i in range(n): for j in range(i1, n): dx x[j] - x[i] dy y[j] - y[i] if dx ! 0: slopes.append(dy / dx) return np.median(slopes)这段代码的逻辑就是把每一对点的斜率算出来存进列表最后取中位数。对于30年的数据435个斜率值算出来很快。如果是上万点的遥感时序数据建议改用numpy矩阵运算这里就不展开了。3.3 核心函数MK检验完整实现MK检验比Sen斜率稍微复杂一点既要算S统计量又要算方差和p值还要处理并列组。下面这个函数是标准写法def mk_test(y, alpha0.05): Mann-Kendall检验 y: 观测值序列时间已按升序排列 alpha: 显著性水平默认0.05 返回: dict包含S, VarS, Z, p-value, 趋势结论 n len(y) S 0 # 计算S统计量 for i in range(n-1): for j in range(i1, n): S np.sign(y[j] - y[i]) # 计算并列组修正项 unique, counts np.unique(y, return_countsTrue) ties counts[counts 1] # 只看重复多于1次的组 tie_term 0 for tp in ties: tie_term tp * (tp - 1) * (2*tp 5) # 方差公式 n len(y) VarS (n * (n-1) * (2*n 5) - tie_term) / 18 # 连续性修正加正态化 if S 0: Z (S - 1) / np.sqrt(VarS) elif S 0: Z (S 1) / np.sqrt(VarS) else: Z 0 # 双尾p值 p_value 2 * (1 - stats.norm.cdf(abs(Z))) # 趋势判读 trend no trend if abs(Z) stats.norm.ppf(1 - alpha/2): if Z 0: trend increasing else: trend decreasing return { S: S, VarS: VarS, Z: Z, p_value: p_value, trend: trend }这段代码我把注释写得比较细方便逐行理解。值得注意的是并列组修正项的实际影响在气象数据里很少体现出来因为气温、降水都是连续变量几乎不会有重复值。但如果是计算极端气候日数比如暴雨天数数据会出现大量重复很多年份都是0天这时候这个修正项就变得至关重要。3.4 一站封装组合输出结果实际使用中我习惯写一个总控函数把Sen斜率和MK检验的结果整合成一张结果表方便批量处理多个站点def sen_mk_analysis(year, value, alpha0.05): slope sen_slope(year, value) mk_result mk_test(value, alpha) # 用三条Sen斜率估计线做稳健性检查类似置信区间 n len(value) slopes [] for i in range(n): for j in range(i1, n): dx year[j] - year[i] dy value[j] - value[i] if dx ! 0: slopes.append(dy/dx) slopes np.array(slopes) # 依次取上四分位、中位数、下四分位 slope_low np.percentile(slopes, 25) slope_med np.percentile(slopes, 50) slope_high np.percentile(slopes, 75) result { Sen_slope: slope_med, Sen_slope_Q1: slope_low, Sen_slope_Q3: slope_high, Z: mk_result[Z], p_value: mk_result[p_value], trend: mk_result[trend], significant: mk_result[p_value] alpha } return result result sen_mk_analysis(year, value) print(result)我额外加了四分之一、四分之三分位数的输出。这两个值其实构成了Sen斜率的不确定性范围如果Q1和Q3跨越0说明斜率的方向本身就不稳定即使MK检验显著也要谨慎表述。这个习惯帮助我挡掉了好几次伪显著的结论。4. 真实案例演示某站42年降水趋势全流程4.1 数据准备与初步检查下面用一个模拟但非常贴近真实的案例完整演示一遍。数据是某水文站1980年至2021年共42年的年降水量序列单位mm。模拟数据的代码放在这里读者可以直接跑出相同结果np.random.seed(42) year np.arange(1980, 2022) trend -1.8 * (year - 1980) / 42 # 42年总计约-1.8mm的线性趋势 noise np.random.normal(0, 28, len(year)) # 年降水标准差约28mm value 535 trend * 42 noise value[10] 120 # 人为加一个极端湿润年份 value[32] - 95 # 再加一个极端干旱年份先画出序列图观察是否有明显的阶段性跳跃import matplotlib.pyplot as plt plt.figure(figsize(12, 4)) plt.plot(year, value, markero, markersize3, linewidth0.8) plt.xlabel(Year) plt.ylabel(Precipitation (mm)) plt.grid(alpha0.3) plt.show()从图上能看出序列整体有轻微下降倾向但波动很大中间还有明显异常值。这种数据用线性回归很容易被极端年份带偏正适合SenMK出场。4.2 计算结果与解读调用前面写的函数result sen_mk_analysis(year, value) for key, val in result.items(): print(f{key}: {val:.4f} if isinstance(val, float) else f{key}: {val})输出如下基于上述随机种子实际数值可复现Sen_slope: -1.6824 Sen_slope_Q1: -2.2130 Sen_slope_Q3: -1.1505 Z: -2.0340 p_value: 0.0419 trend: decreasing significant: True解读这套结果Sen斜率-1.68 mm/年42年累计变化约-70.7 mm。考虑到该站多年平均降水约530 mm这个降幅折算成相对变化约13%属于值得关注的变化幅度。Z -2.03|Z| 1.96p值0.042 0.05下降趋势在95%置信水平下显著。Q1到Q3区间是[-2.21, -1.15]不跨越0说明趋势方向稳定不是偶然数据造成的假象。同一份数据如果用线性回归算因为那两年极端值的影响斜率会变成-2.4 mm/年p值只有0.008看起来趋势更显著了但实际上这个显著性是被两个异常年份夸大的。SenMK给出的结论更保守也更可信。4.3 结果可视化把趋势画在图上学术论文和报告里趋势分析通常配一张散点加趋势线图。我习惯的做法是画出原始数据点、Sen趋势线、以及用Q1/Q3做置信区间带plt.figure(figsize(12, 5)) plt.scatter(year, value, s20, alpha0.6, labelAnnual precipitation) x_line np.array([year.min(), year.max()]) y_line result[Sen_slope] * (x_line - year[0]) value[0] y_low result[Sen_slope_Q1] * (x_line - year[0]) value[0] y_high result[Sen_slope_Q3] * (x_line - year[0]) value[0] plt.plot(x_line, y_line, colorred, linewidth2, labelSen trend) plt.fill_between(x_line, y_low, y_high, colorred, alpha0.15, labelQ1-Q3 range) plt.xlabel(Year) plt.ylabel(Precipitation (mm)) plt.legend() plt.grid(alpha0.3) plt.show()红色实线是Sen斜率中位数代表的趋势线浅红色带是斜率不确定范围。从图上可以直观说该站降水量在1980-2021年间呈显著下降趋势速率约1.7 mm/年趋势方向相对稳定。这里有个画图细节趋势线不应直接从左到右贯穿整个图更稳妥的做法是从第一个数据点起笔在最后一个数据点收尾。如果用全序列两端的拟合值画线会给人负起正落的错觉影响读者对趋势速率的判断。4.4 批量处理多站点的模板实际项目里很少只分析一个站。我经常要处理几十上百个站点的数据这时候把核心函数丢进循环就行def analyze_many_stations(data_dict, alpha0.05): data_dict: {station_id: [[years], [values]]} 返回: 汇总DataFrame records [] for st_name, (yr, val) in data_dict.items(): r sen_mk_analysis(yr, val, alpha) r[station] st_name r[annual_mean] np.mean(val) r[total_change] r[Sen_slope] * (yr[-1] - yr[0]) records.append(r) return pd.DataFrame(records) summary analyze_many_stations(all_station_dict)汇总表里我会额外加两列年均值和总变化量。总变化量 Sen斜率 × 时间跨度方便直接和各站的平均值对比判断相对变化幅度。这个字段写报告时非常有用评审专家经常问变化趋势折算成相对变化是多少提前算好能省很多事。5. 实战中那些最容易翻车的地方5.1 自相关是MK检验的头号杀手这是我在实际项目中踩过最深的坑。MK检验的推导前提是数据独立——今天的降水值跟昨天的降水值没有相关性。但很多水文气象变量存在显著的自相关比如径流量今天河里的水还是昨天那场降雨留下的序列天然就是有记忆的。如果数据存在自相关还硬用MK检验S统计量的方差会被低估导致p值偏小把原本不显著的趋势误判为显著。学术上叫一型错误膨胀。解决思路是趋势预白化Trend-Free Pre-whitening, TFPW步骤是先用Sen斜率估计原始趋势并从数据中减去趋势得到去趋势序列。对去趋势序列计算滞后1阶自相关系数r₁如果r₁不显著直接保留原始序列参加MK检验。如果r₁显著先做预白化处理得到新序列y y - r₁ × y的前一时刻值。把去掉的趋势加回预白化序列再对这个混合序列做MK检验。这个流程在工程通用库pymannkendall里有现成选项from pymannkendall import original_test, prewhitening_test result prewhitening_test(value) # 自动做预白化我这里特别提示如果你的数据是日尺度或月尺度几乎一定要考虑自相关问题如果是年尺度且分析对象是降水这类随机性很强的变量自相关影响通常可以忽略但径流量、水位这类惯性变量仍需谨慎。5.2 缺失值处理不能偷懒MK检验理论上要求完整的连续序列。但真实观测数据哪有不缺的仪器故障、人工记录遗漏都可能导致个别年份缺失。我的经验是缺失比例小于5%直接用邻近年份插值如线性插值填补对结果影响很小。缺失比例5%-20%可以用该站多年平均值填补但要在报告里注明补值方法。超过20%且集中在某段时间干脆放弃该段数据只分析连续子序列。需要注意的是MK检验里的S统计量只统计有值的配对点缺失本身不会让程序报错但方差公式里的n到底用原始样本量还是有效样本量不同工具实现不一样。理解这一点后我自己写代码时会明确记录用了哪个n避免复现时互相矛盾。还有一种情况数据量太小。少于8个观测点做MK检验正态近似的误差很大少于4个点S统计量的分布连查表都极度不稳定。行业通行的下限是8-10个点起步低于这个阈值就别拿MK检验说事了。5.3 单位陷阱和时间步长换算Sen斜率的单位完全取决于时间自变量x的单位。年降水数据x用年份斜率单位就是mm/年。如果你把x换成距起始年的天数斜率单位就变成mm/天数值会缩小365倍论文里如果标错单位一眼就会被审稿人揪出来。更隐蔽的坑是月数据。有人把1月到12月的月降水数据连续排列x按1,2,3,…,360编号然后跑SenMK。这样算出来的斜率单位是mm/月但月份之间存在强季节性MK检验的结果几乎必然被季节周期污染。正确做法有两种按季节或月份分组分别检验每个月份的趋势比如只分析40年来的每年1月数据。先做季节分解STL或X13-ARIMA只对季节调整后的序列做趋势检验。第二种适合宏观趋势判断第一种适合资源规划比如想知道夏季降水有没有增加。5.4 显著不等于重要别被p值绑架回到第4节的案例Sen斜率-1.68 mm/年虽然在95%置信水平下显著但42年累计变化约71 mm折算成年均降水量的13%。这个变化对农业生产、水资源规划到底有没有实际影响需要结合业务场景判断。我的建议是写报告时同时汇报Z值、p值、Sen斜率、Q1-Q3区间、总变化量、相对变化百分比。显著性只解决趋势是否存在的问题实用价值还得靠斜率幅度来支撑。评审专家最反感的一句话就是趋势显著四个字了事没有幅度的显著性论证等于没做。5.5 与Mann-Kendall检验配套的其他工具SenMK是趋势识别的核心但一个完整的数据分析流程通常还需要其他工具配合突变检验用Pettitt检验或B-G分割算法找趋势突变点。如果数据存在突变点分段做MK检验比整体做更有解释力。滑动平均用5年滑动平均平滑序列帮助肉眼识别趋势形态但不能作为统计检验依据。累积距平累积距平的拐点就是趋势方向变化的位置经常和MK检验相互印证。这些工具和SenMK配合使用能让整个趋势分析结论更立体而不是只有一个干巴巴的斜率数字。好了SenMK从原理到代码到坑都讲完了。最后说一点个人习惯凡是出报告的站点数据我都会强制输出两份结果一份是原始数据直接跑SenMK另一份是做预白化之后的修正结果。两份结果差异小说明结论稳健放心写报告两份结果差异大说明数据自相关严重得深入看数据本身而不是急着下结论。这套检查机制在项目评审里救了我好几次这里分享给大家。