Sen+MK趋势分析实战：Python实现水文气象时间序列稳健趋势检验

1. 当你不确定数据到底在涨还是在跌的时候就该拿出SenMK了做水文、气象、环境监测这类长时间序列数据分析的人迟早会遇到同一个问题手里拿着一二十年的观测数据曲线图出来以后肉眼看着好像有点变化但你拿给领导看、写进报告里、或者发论文的时候不能只说“我觉得它在上升”。你得拿出统计证据告诉别人这个趋势是真的存在而且幅度到底有多大。我以前处理泾河流域某个水文站的年径流数据时就踩过这个坑。直接用最小二乘法拟合线性趋势R平方只有0.12p值倒是有显著性但画出来残差分布是歪的数据里还有几个极端枯水年的离群点抓着回归线的尾巴不放。后来换了Sen斜率估计配合Mann-Kendall检验结论一下子清爽了趋势方向、显著性、变化速率全都清清楚楚而且对那几个离群点几乎免疫。这套组合在学界和工程界都叫“SenMK趋势分析”本质上就是两种非参数方法的搭配使用Sen斜率Theil-Sen Estimator负责估计趋势的强度也就是每单位时间变化多少Mann-Kendall检验负责判断趋势有没有统计显著性也就是这个趋势到底是真信号还是随机波动。核心优势是它们都不要求数据服从正态分布也不怕少数异常值的干扰对缺失值和非线性趋势也有不错的容忍度所以成了趋势分析领域最常用的组合之一。这篇文章不打算堆公式推导就从一个实际使用者的角度把SenMK的原理、Python实现、结果解读办法、以及我在实际项目中踩过的问题一次说清楚。如果你手头正好有一批环境监测数据、水文气象数据或者长时序的业务指标数据这套方法大概率能直接帮你出结论。适合谁来参考一类是刚开始接触趋势分析的研究生和科研助理需要一种稳健、有据可查的分析工具另一类是环保、水利、农业领域的从业人员需要在报告里给出趋势判断和量化速率。文章里我会给出可复现的Python代码和完整步骤你拿去改改路径就能跑。2. 为什么是它们俩非参数方法的底气从哪来2.1 最小二乘法拟合趋势为什么不总是好用先说传统思路的问题否则你不理解SenMK的价值。常规做法是对时间t和观测值y做线性回归得到斜率b然后用回归系数的p值判断趋势是否显著。这个套路在数据干净、服从正态分布、方差齐性的情况下确实没问题但真实的环境监测数据哪有这么规矩第一个问题是离群值。某个年份赶上极端天气数据一下偏离主群三四倍标准差最小二乘的拟合线会被强行拽过去导致斜率估计严重失真。第二个问题是数据分布。降雨、径流、污染物的浓度数据往往是有偏分布有的是正偏态有的是对数正态回归模型里的正态性假设经常不成立。第三个问题是缺失值。野外监测站不可能保证每个时间段都有完整记录仪器故障、工况调整、采样计划变化都会造成缺口线性回归对样本空白的处理比较生硬。我见过不止一个案例用最小二乘法估计出来的趋势方向是正的p值也没问题但把异常年份去掉之后再做一遍结论直接反转。这说明你对离群点太敏感了稳健性不行。SenMK恰恰在这些场景下表现稳定因为它们不依赖数据的分布形态天然抗少数离群值的干扰。2.2 Sen斜率所有点对斜率的中位数天然抗噪声Sen斜率估计的思想非常直观一句话就能说明白把数据里任意两个点之间的斜率都算一遍最后取所有斜率的中位数这个中位数就是整个序列的趋势速率。为什么取中位数而不是平均值因为中位数是分位数的一种天然对极端值不敏感。假如你有10000个点对斜率里面有几十个因为离群点产生了极端值取平均容易把结果带偏但取中位数就稳稳地落在中间位置几乎不受影响。这就是稳健统计里最常见的思路用秩代替具体值用分位数代替均值。具体形式上对于时间序列 \((t_i, y_i)\)任取两个不同的时间点 \(i\) 和 \(j\)计算斜率\(slope_{ij} \frac{y_j - y_i}{t_j - t_i}\)对所有 \(i j\) 的组合都算一遍最后取这些斜率的中位数作为Sen斜率记为 \(\beta\)。这个 \(\beta\) 的含义就是“单位时间内目标变量平均变化了多少”比如每年增加多少毫米的降水量、每年升高多少摄氏度的气温、每年下降多少立方米的径流量。单位由你的输入数据决定非常直观。有人会问如果时间序列很长点对数量会不会爆炸比如n365的逐日序列点对数量就是 \(n(n-1)/2 \approx 66000\) 个计算量确实不小但现代计算机完全扛得住。如果你的数据是逐分钟的n上十万就要考虑向量化实现或者抽样计算了后面代码部分我会给出高效写法。2.3 Mann-Kendall检验用秩次变化判断趋势是否显著Sen斜率告诉你趋势有多快但没告诉你这个趋势是不是“真”的。如果一个序列本身就是随机游走没有真实的单调趋势你用Sen斜率去算它照样能给一个数值出来。所以必须加一道显著性检验的关卡这就是Mann-Kendall检验的活儿。Mann-Kendall检验的基本思想是把序列里所有前、后的数据对拿出来比较统计有多少对是“后一个比前一个大”上升趋势多少对是“后一个比前一个下降趋势”的。如果上升的数目和下降的数目差不多说明数据随机波动没有系统性趋势如果某一方向的数量明显压倒另一方向说明序列里存在显著的趋势信号。算出来的核心统计量叫S统计量再标准化成Z值最后用标准正态分布算出p值。p值小于0.05我们一般认为趋势显著小于0.01则称极显著。这里需要留意一点MK检验的原假设是“序列不存在单调趋势”备择假设是“存在单调趋势”所以p值越小越有理由拒绝原假设也就可以认为趋势是真实存在的。这两个方法天然互补Sen斜率负责告诉你是升是降、速率多少MK检验负责告诉你这个升降能不能从统计上站住脚。在实际汇报或者论文中标准的写法就是“MK检验表明该序列在α0.05水平上存在显著上升趋势Sen斜率为每年0.32毫米。”简洁、准确、有据可查。2.4 为什么这套组合被称为“趋势分析的黄金搭档”单独用其中某一个其实都不够完整。只用MK检验你只能知道“有趋势”或者“没趋势”但没法告诉别人趋势有多强写报告时少了关键的定量指标只用Sen斜率你能算出速率但万一这个速率是从一堆随机波动里算出来的你会被质疑结论的可信度。所以行业里约定俗成地把两者打包使用。MK检验给结论Sen斜率给参数两者一起构成完整的趋势分析输出。任何一个专业的水文气象分析软件里比如美国的趋势分析工具、水利系统的相关规范都能看到这个组合的影子。你写论文的时候审稿人也默认接受这套方法的有效性因为它不依赖具体分布假设适用范围广解释起来也没有太大的争议。3. 手把手实操从原始时间序列到趋势结论的完整流程3.1 数据准备先搞清楚什么格式才能喂给算法SenMK的输入非常朴素就是一组按时序排列的观测值。你可以准备两列数据一列是时间年份、月份、日期都可以一列是对应的观测值。但要注意一点时间步长最好保持一致要么都是年值要么都是月值尽量不要把不同频率的数据混在一起分析否则取点对斜率的时候不同位置的时间间距不一致会干扰结果解释。我一般建议分析之前做三步预处理。第一步是缺失值检查如果缺失比例小于5%可以用线性插值补齐如果缺失比例超过10%或者存在连续的大段空白就要慎重了因为MK检验的方差公式里默认序列是连续单调的大段空白会影响方差估计的准确性。第二步是异常值筛查这里说的异常值不是指极端年份的真实记录而是指明显的记录错误比如流量出现负值、降水量超过物理上限这些要先订正。第三步是季节性处理如果你的数据是月尺度或者日尺度并且存在明显的季节周期性直接套用原始MK检验会把季节信号误判为年际趋势需要改用季节版MK检验或先做季节分解这个我在第5部分详细说。3.2 手写实现核心计算搞懂内部逻辑再上库函数我之前带过不少学生和同事发现一个规律直接用现成的库一旦结果看起来不太对很多人就完全不知道从哪里排查了。所以我建议至少先手写一遍核心算法自己搞清楚每个数字是怎么来的然后再切换到成熟的库。下面是我写的一个简洁版Sen斜率估计函数思路清晰适合用来理解原理import numpy as np def sen_slope(t, y): 计算Sen斜率Theil-Sen斜率估计 参数 ---------- t : array_like 时间序列的时间变量年份/月份等 y : array_like 时间序列的观测值 返回 ------- float Sen斜率的中位数估计 t np.asarray(t) y np.asarray(y) n len(y) slopes [] for i in range(n - 1): dt t[i1:] - t[i] dy y[i1:] - y[i] # 避免除零 valid dt ! 0 slopes.extend(dy[valid] / dt[valid]) return np.median(slopes)这个版本使用了NumPy的向量化计算避免了双重循环带来的性能开销对几千个数据点完全够用。如果你有十万级的数据量还可以用scipy.stats.theilslopes直接调用现成实现内部逻辑是一样的但它还会给你返回置信区间这个很有用。然后是Mann-Kendall检验的手写版本from scipy.stats import norm def mann_kendall(y): Mann-Kendall趋势检验不考虑季节性与自相关 参数 ---------- y : array_like 时间序列观测值 返回 ------- dict 包含S统计量、方差、Z值和p值的字典 y np.asarray(y) n len(y) # 计算S统计量 s 0 for i in range(n - 1): diff y[i1:] - y[i] s np.sum(np.sign(diff)) # 处理并列值ties对方差的影响 unique, counts np.unique(y, return_countsTrue) ties np.sum(counts[counts 1] * (counts[counts 1] - 1) * (2 * counts[counts 1] 5)) # 计算方差 n_f float(n) var_s (n_f * (n_f - 1) * (2 * n_f 5) - ties) / 18 # 计算标准化Z值带连续性修正 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值 p 2 * (1 - norm.cdf(abs(z))) return {S: s, variance: var_s, Z: z, p_value: p}代码里有一个细节值得解释ties的处理。当数据中出现重复值的时候S统计量里那些等于0的差值不贡献方向信息但会减少方差如果不在方差公式里做修正p值会偏低容易把不显著的趋势误判为显著。所以只要数据有并列值必须做这个修正别偷懒。3.3 更省事的现成方案用pymannkendall一行出结果自己手写很好但真正跑项目的时候直接用成熟的库更稳妥因为库作者对各种细节的处理比我们临时写一遍想得更周全。Python生态里处理MK检验最方便的是pymannkendall这个库它支持原始MK检验、季节MK检验、相关趋势检验等多种变体接口非常统一。安装和基本用法如下pip install pymannkendallimport pymannkendall as mk import numpy as np # 假设 data 是一维数组按时间顺序排列 data np.array([12.3, 12.8, 11.9, 13.1, 13.5, 13.2, 14.0, 13.8, 14.5, 14.9]) # 原始MK检验 result mk.original_test(data) # 输出主要结果 print(趋势方向:, result.trend) print(p值:, result.p) print(Sen斜率:, result.slope) print(截距:, result.intercept) print(Z值:, result.z) print(S统计量:, result.s)这个库的输出结果里有几个关键字段需要解释一下。trend是趋势方向可能是increasing、decreasing或no trendp是显著性水平slope是Sen斜率单位跟随你的输入数据如果输入的是逐年数据单位就是“每单位时间”的变化量intercept是趋势线的截距配合斜率可以用来画趋势线。实际工作中我一般不会直接跑一遍原始MK就结束通常会先画图看数据形态再用多条MK变体交叉验证。pymannkendall里还提供了seasonal_test和correlated_series_test等函数对处理季节性和自相关数据很有帮助。3.4 结果解读别只盯着p值斜率置信区间同样重要跑完代码拿到结果很多新手会犯一个毛病只看p值小于0.05就说“趋势显著”其他输出一概不管。这样做是有隐患的。p值告诉你统计上是否显著但显著不等于有实际意义尤其是在样本量很大的时候一个0.001的微小变化也可能被标成显著但它在业务上可能毫无价值。我的建议是看三个东西。第一是趋势方向这个最直观是上升还是下降先确认跟你实际观测的直觉是否一致如果不一致先找数据问题不要急着出结论。第二是Sen斜率的大小结合业务背景判断单位时间的变化量有没有实际意义比如径流每年减少0.01立方米每秒统计上可能显著但对水库调度几乎没有影响。第三是斜率的置信区间scipy.stats.theilslopes这个函数会返回low_slope和high_slope也就是斜率的95%置信区间区间越窄说明斜率估计越稳定区间跨越0说明趋势其实在统计上不显著跟MK检验的p值互相印证。还有一点需要提醒MK检验的单数检验和双侧检验需要提前搞清楚。上面代码算出来的是双侧检验的p值默认用来判断“是否存在任意方向的趋势”。如果你有明确的方向性假设比如只关注下降趋势那要用单侧检验的p值。这些细节在论文方法部分写清楚不然审稿人较真的时候你会很被动。4. 真实案例某水文站年径流序列的SenMK分析全记录4.1 案例背景与数据说明这套分析从原理到代码都讲完了接下来用一个完整的案例把整个流程串起来。我的数据取自某水文站1990年到2020年的年径流量记录一共31个样本点单位是亿立方米。这类数据有几个典型特征年际波动大偶尔出现特枯年和丰水年数据分布通常右偏而且可能存在一阶自相关。这些特征决定了它非常适合用SenMK组合来分析。拿到数据的第一件事永远是画图。把年径流量按时间轴画成散点图这里不加任何趋势线只是观察形态看有没有明显的趋势方向、突变点、周期波动。4.2 一步一步走完整套流程第一步做描述性统计分析看数据的集中趋势和离散程度import pandas as pd import numpy as np from scipy import stats # 构造数据示例实际使用时替换为你的数据 years np.arange(1990, 2021) q np.array([ 96.2, 102.5, 88.7, 110.3, 98.4, 92.1, 105.8, 99.6, 87.3, 115.4, 103.2, 91.0, 95.8, 99.1, 87.6, 93.5, 86.2, 101.4, 97.3, 84.5, 108.6, 95.4, 90.2, 105.7, 94.8, 88.2, 102.9, 93.1, 89.5, 107.2, 92.8 ]) print(f样本数量: {len(q)}) print(f均值: {np.mean(q):.2f} 亿立方米) print(f标准差: {np.std(q):.2f} 亿立方米) print(f变异系数: {np.std(q) / np.mean(q):.3f})第二步计算Sen斜率和截距我这里同时用我自己手写的函数和scipy.stats.theilslopes交叉验证from scipy.stats import theilslopes # 自定义Sen斜率 beta_custom sen_slope(years, q) print(f自定义Sen斜率: {beta_custom:.4f}) # scipy实现顺便获取置信区间 res theilslopes(q, years, 0.95) print(fSciPy Sen斜率: {res.slope:.4f}) print(f95%置信区间: [{res.low_slope:.4f}, {res.high_slope:.4f}]) print(f截距: {res.intercept:.4f})第三步执行MK检验# 手写版本 mk_result mann_kendall(q) print(f手写MK: S{mk_result[S]}, Z{mk_result[Z]:.3f}, p{mk_result[p_value]:.4f}) # pymannkendall库版本 import pymannkendall as mk_lib mk_lib_result mk_lib.original_test(q) print(f库MK: 趋势{mk_lib_result.trend}, Z{mk_lib_result.z:.3f}, p{mk_lib_result.p:.4f})第四步把趋势线叠加到原始时间序列图上直观呈现变化情况。这里用的趋势线就是Sen斜率 × 时间 截距画出来一条直线方便可视化对比。4.3 结果怎么解读才够专业像我构造的这份示例数据运行完以后结果大概是这样的Sen斜率是负值大约-0.15意味着年平均径流量以每年0.15亿立方米的速度减少31年间总共减少约4.6亿立方米这个量级占均值的5%左右只能算轻度下降。MK检验的p值在0.2附近大于常用的0.05显著性水平说明这个下降趋势在统计上并不显著也就是说它有可能只是长期气候波动中的一个低频振荡。这时候报告里不能写“径流量显著下降”而要写成“1990—2020年径流量呈不显著下降趋势p0.05Sen斜率为-0.15亿立方米/年”。这个表述虽然听起来没那么“厉害”但从统计角度是严谨的经得起推敲。科研和工程报告中最常见的错误就是只看趋势线下降就像下结论说“减少显著”实际上连MK检验都没跑过。反过来也有一些场景是趋势线几乎水平但MK检验p值反而很小这说明存在某种非线性但系统性的一致方向性变化线性回归看不出来MK检验能捕捉到。所以每一步都有它不可替代的价值。5. 实战避坑指南这些坑我替你踩过了5.1 有季节性周期的数据直接用原始MK会出错我之前接过一个项目分析某个断面逐月的水质监测数据用原始MK检验跑出来“显著上升趋势”结果画月度时间序列图一看水的就是季节性周期在作怪夏天浓度高、冬天浓度低MK检验把每个季节内部的升降全都算了一遍最后合并出了一个整体向上的假趋势。不是水质真的在变差而是数据的季节性没有被剔除。解决办法有两个路径。第一个路径是在分析前做月份逐月分离比如1月和1月比较、2月和2月比较然后用季节版MK检验seasonal Mann-Kendall来处理。pymannkendall里直接有seasonal_test函数。第二个路径是先做季节分解或者去季节性处理比如用移动平均或STL分解提取出趋势分量然后对趋势分量做常规MK检验。两条路都可行我个人的习惯是优先用季节版MK因为它的原假设和备择假设设计得更严密而且操作上不容易引入二次处理的偏差。年际序列比如一年一个值的年径流、年降水总量一般不需要考虑这个问题但月尺度和日尺度的数据几乎必然要处理季节性这点切记。5.2 自相关数据会让p值“假显著”注意先做预白化时间序列数据最烦人的特征除了季节性是自相关。比如逐日气温数据今天30度明天大概率也在30度附近这种数据不满足MK检验里“样本独立”的假设。自相关会让MK检验的方差被低估进而导致p值被人为压低把没有显著趋势的序列误判成有趋势也就是“假显著”问题。如果你的数据存在明显的一阶自相关步骤上应当先做预白化处理也就是用AR(1)模型把数据的自相关结构去掉再对残差做MK检验。pymannkendall里对应的函数是test_inequality和correlated_series_test后者专门处理具有自相关结构的时间序列。如果数据有更高阶的自相关建议先拟合AR(p)模型然后对残差做趋势分析。如何判断数据有没有自相关最直接的办法是算一下序列的Lag-1自相关系数或者画ACF图。如果Lag-1自相关系数的绝对值超过0.2就要小心了。这个指标不高但很多人忽略它实际分析中会因为自相关而翻车的情况远比想象的要多。5.3 样本量太小MK检验的近似正态分布会失真MK检验在样本量较大的时候用Z值近似标准正态分布来计算p值这个近似在n30的时候效果不错。但如果只有几年数据比如n10或更少正态近似的误差会变大置信度就不太够。面对小样本有两种处理思路。第一种是走精确检验路线直接枚举所有可能的排列组合计算精确的p值基础代码里可以自己实现或者使用一些统计软件里的非参数检验功能。第二种是老老实实承认结论的局限性在报告里注明“样本量有限检验效能较低”。我自己做项目时如果n小于20我会把MK结果当作一种参考信号而不是定论需要更多年份的数据积累后再做最终判断。这不是回避问题做数据分析本来就要对结论的边界条件诚实。5.4 多个站点多个序列同时分析时记得做多重检验校正这个坑做单序列分析时碰不到一旦你同时分析十个站点的趋势每个站点独立做MK检验就要考虑多重检验的问题。假设显著性水平是0.05纯随机情况下十个站点里也会有约0.5个站点被误判为“显著”多站点分析时误判概率会被放大。我处理过一个流域内20个子流域的降水趋势独立跑MK检验出来3个站点显著看着好像有空间规律但一做Benjamini-Hochberg校正剩下了不到1个站点显著。这是很现实的差距。如果你做的是多站点对比研究建议对p值做FDR错误发现率校正做法很简单用Python里的statsmodels.stats.multitest.multipletests对p值列表做一次FDR校正即可。5.5 可视化的细节趋势线别乱画显著性也要标清楚最后说一个偏颜值和表达的话题。趋势分析结果画图的时候很多新手会把原始数据点和Sen趋势线一起画出来却在图上不标注显著性。读者看到一条下坡线就以为是显著下降但你明明p值都0.3了这条线根本不具备统计显著性。所以图上的建议标注方式是在标题或者图例里明确写“趋势不显著”或“趋势显著p0.05”让看图的人不会被误导。画趋势线的时候不要用最小二乘回归线的斜率来画要用Sen斜率。很多可视化工具默认做线性拟合你得手动指定斜率和截距。用Matplotlib画的时候可以直接用Theil-Sen估计出来的斜率乘以时间加截距来画直线代码层面就是import matplotlib.pyplot as plt slope res.slope intercept res.intercept trend_line slope * years intercept plt.figure(figsize(10, 5)) plt.scatter(years, q, s40, label观测值) plt.plot(years, trend_line, r--, labelfSen趋势线 (斜率{slope:.2f})) plt.xlabel(年份) plt.ylabel(径流量亿立方米) plt.title(年径流量趋势分析) plt.legend() plt.grid(alpha0.3) plt.show()这里的标题可以先加一个占位符等MK检验结果出来之后手动把显著性信息补上。图是报告的脸面数据再漂亮图乱七八糟也白搭。最后几件小事SenMK这套组合我用过无数次从科研论文到工程报告都跑过。它最大的优点是皮实不管数据来源是气象站、水文站还是环境监测站不管数据量多少只要按步骤来基本都能得到可信的结论。但它也不是万能的它假设趋势是近似单调的如果你的序列存在明显的突变点、周期性大起大落或者非线性变化比如U型变化MK检验可能会把它们识别为“无显著趋势”这时候建议配合突变检验比如Pettitt检验、Mann-Kendall突变点检验一起使用。分析趋势这件事方法和结论之间永远隔着你对数据的理解。任何统计量都替代不了对业务背景的把握——在解释结果前先搞清楚你的数据是怎么来的、经历了哪些变更、代表什么物理意义。这样跑出来的趋势结论才是真正能说服人的结论。如果你在实际分析中遇到了奇怪的结果多半不是数学出了问题而是数据在某个环节埋了坑。这时候不妨回到原始记录重新捋一遍。这套方法的容错能力很强但前提是你用对了地方、看全了输出。