搞地学、遥感、水文、气象数据分析的老哥老姐们SenMK趋势分析这组词估计十有八九都熟。它基本是“时间序列趋势检测”里的默认组合了Sen负责估算变化速率MK负责判断趋势显著性两个一配合既能告诉你“涨了还是跌了”又能告诉你“这个涨跌能不能信”。我自己日常处理的不是NDVI就是降水、土壤湿度凡是论文里写“分析近二十年演变趋势”的基本都用这套方法交差。这篇文章就把我从原理理解、数据准备到Python跑批的完整流程连同踩过的坑一起拆开给刚开始接触趋势分析的学生也给每一个被长序列数据折磨的工程师。先给个结论这套方法不挑数据分布、对离群值不敏感、实现成本低属于那种“你不需要懂太多数学也能用对”的工具。但它也有边界比如对自相关数据、季节数据、短序列数据直接套用会出大问题。下面我会先把原理讲透再给完整代码和参数选择思路最后把高频问题列成排查清单照着做基本能避掉九成以上的坑。1. 为什么是“SenMK”这对组合解决什么问题1.1 常规线性回归的尴尬被离群值牵着走拿到一条时间序列多数人第一反应是做一元线性回归看斜率是否显著。这在教科书上没问题真实数据却常常不配合某一年降水特别大、某一年NDVI因为传感器问题突然跳变、或者源数据本身就带着几个极端噪声点。最小二乘回归对这类离群值非常敏感一个异常点就足以把斜率从正拽成负或者从显著拉成不显著。我印象很深的一次教训某站点20年降水序列我做线性回归得到的趋势是“每年增加约3.1毫米p0.02”看上去很漂亮。后来回查原始记录发现第7年有一次极端特大暴雨数值比正常年份高出好几倍。把这个点去掉再算回归斜率直接变成负的。一个点翻转整个结论这显然不合理。后来又用Sen斜率一算数字稳定在“每年下降约0.8毫米”因为Sen拿的是所有点对斜率的中位数那个极端暴雨点最多影响它参与的少数几个点对不足以撼动中位数。这也是为什么在真实项目中我宁可把线性回归当作初筛最终结论一律以SenMK为准。如果你不理解中位数的意义可以这样类比一群人里混进一个“年薪百亿”的算平均收入瞬间被拉爆但算收入中位数大家的生活水平几乎不受影响。Sen斜率就是“收入中位数版”的趋势估计它不追求完美拟合只求给一个不被异常值绑架的稳健方向。1.2 Sen斜率估计把所有点对斜率排个队Sen斜率也叫Theil-Sen估计计算起来非常朴素对于一条长度为n的序列任意取两个时间点i和ji 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分所有分数加总得到统计量SS Σ 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)(2n5) - Σ_t t_i(t_i-1)(2t_i5)] / 18后一项是“并列值修正”。如果数据里有大量完全相等的点比如降水经常出现0.0毫米的月份这个修正必须做否则方差被低估p值虚小容易得出“假显著”结论。Z值近似服从标准正态分布所以双侧检验p值就是2倍的尾部概率p 2×(1 - Φ(|Z|))。p小于0.05我们就说趋势显著小于0.01说极显著。MK检验不要求数据正态不怕离群值对大多数单调趋势的检测都很有效所以成了环境变化研究里的常青树。1.4 为什么两者总是成对出现Sen和MK能经常合体是因为它们优势互补、方法论基因一致。工具回答的问题输出Sen斜率变化速率是多快方向是增还是减斜率β单位是“每年多少量”MK检验这个趋势是否统计显著Z值、p值、Kendall TauSen算出“每年增加0.003的NDVI”MK负责告诉你“p0.42这个增加很可能是噪声”。反过来MK说“z2.31p0.02存在显著趋势”Sen则告诉你是每年1毫米还是每年50毫米。只用一个都容易误读一个“显著”但没有量级一个“有量级”但不可信。在实际项目中我还会顺带输出MK计算过程中得到的Kendall Tau它介于-1和1之间表示趋势的强度。写报告时表格里放“Sen斜率、p值、Tau值”三列基本就齐活了。2. 动手前的准备数据形态、工具选型与几项关键约定2.1 数据应该是什么样一维时间序列和多波段栅格SenMK的输入按项目场景分两种。第一种是单点时间序列比如某个气象站的逐年降水量或者某个采样点的水质指标一维数组长度为年份数。第二种是空间栅格比如MODIS NDVI的年合成产品20年就是20个波段每个波段代表一年的空间分布合成一个多波段GeoTIFF。处理栅格时本质上是把每个像素对应的“时间维数组”单独拎出来做一次SenMK再把结果写回成一个单波段栅格。无论哪种形态有一个前提必须搞清楚时间分辨率是什么。年度数据用普通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 matplotlibpymannkendall这个库非常完整原始MK、季节MK、自相关修正MK、预白化MK全都有函数名也直观。如果你不用PythonR语言里也有zyp和trend包功能类似但代码风格不同。这里我只讲Python路径毕竟栅格批处理这一块Python生态确实更顺手。2.3 两个常用实现路径手搓公式与现成函数做趋势分析最忌讳的就是黑盒使用。我的习惯是先用现成函数快速出结果再手写一遍核心公式验证理解。项目紧的时候直接调库项目松的时候至少把S统计量和Sen斜率自己算一遍这样出了意外结果才知道往哪个方向排查。pymannkendall最核心的接口是original_test它返回一个命名元组里面包含trendincreasing/decreasing/no trend、h是否拒绝原假设、pp值、zZ统计量、TauKendall Tau、sS统计量、var_sS的方差、slopeSen斜率、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_valid15): series np.asarray(series, dtypefloat) 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, dtypefloat) 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)你可能会担心这个双重循环太慢。确实n20时有190个点对n50时则上升到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(fSen斜率: {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, dtypefloat) 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(fS {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 Tauoriginal_test返回的slope就是第3.2节手写函数的结果两者应该对得上。注意result.trend对“是否拒绝原假设”的判断使用的是0.05显著水平如果你想要0.1或者0.01别直接用trend字段改用p值自己做判断。3.4 第四步批量处理栅格逐像素确实慢但别怕真正让SenMK项目从“练习”变成“实战”的是空间栅格批处理。假设你有一个NDVI年合成多波段GeoTIFFshape是(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(count1, dtypefloat32, nodatanp.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(fprocessed {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值最终画图时还要再分类。我常用的分级规则如下分级判定条件含义极显著上升Sen 0 且 p 0.01趋势非常可靠显著上升Sen 0 且 0.01 ≤ p 0.05趋势可靠弱显著上升Sen 0 且 0.05 ≤ p 0.1仅作参考趋势不显著p ≥ 0.1不能判断有趋势弱显著下降Sen 0 且 0.05 ≤ p 0.1仅作参考显著下降Sen 0 且 0.01 ≤ p 0.05趋势可靠极显著下降Sen 0 且 p 0.01趋势非常可靠出图时我习惯用暖色表示上升、冷色表示下降不显著区域统一用浅灰或透明。如果你处理的指标有物理范围比如NDVI理论范围是-1到1别忘了在色带数值上标注清楚单位比如“Sen斜率: 0.001/年”避免读者误读成“每年增加1%”。另外我还会生成一张“有效样本量”栅格当作辅助诊断图。当某些区域趋势看上去很奇怪时第一件事就是查那里的有效年份数是不是太少。很多时候一个显著性“突变区”其实是缺测造成的假象有了样本量栅格这种问题一目了然。4. 问题排查与避坑实录4.1 有效样本不足和NaN趋势结果不可信这是栅格处理中最常见的坑。MODIS NDVI常因为云覆盖导致某些年份某些像素为无效值降水数据遇到传感器故障也可能缺测。如果缺失年份太多MK检验的方差估计会失真。举例来说20年数据里只有8年有效序列长度和“自由度”都不足即使真的存在趋势也很难检测出来。我的做法是给每个像素设一个阈值比如有效样本≥15才返回结果低于阈值直接输出NaN并mask掉。这个阈值不是拍脑袋而是从检验功效角度考虑的n15时MK对中等强度的趋势尚能检测n8时基本什么都检不出来。如果你的时间序列本身就是短序列比如只有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里带自相关修正的方法。这个库提供了几个现成接口方法适用场景original_test独立同分布数据常规年度序列hamed_rao_modification_test方差修正版适合存在自相关的数据yue_wang_modification_test另一种方差修正方式trend_free_prewhitening_modification_test先去趋势再预白化适合正自相关明显的数据我的经验阈值是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)同理季度数据用period4周数据用period52。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(count1, dtypefloat32, nodatanp.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(windowwindow) # shape: (n_time, block_h, block_w) # 对block做像素级SenMK再把结果写回对应窗口计算过慢则取决于你点对循环的效率。如果年份数超过40每个像素点对数量至少780对全国范围内就是几十亿次比较。我的建议是核心函数用numba的njit装饰或者用numpy的向量化矩阵构建点对差分以牺牲部分内存换速度再或者用multiprocessing多进程并行。对单个项目而言先用几万像素小范围测试正确性再上全图避免程序跑了两小时才发现结果全是NaN。4.5 p值显示和解读别掉进“科学计数法”的坑MK输出的p值有时候是1e-16这种极小数字很多人偷懒写成“p 0.000”。这个习惯在正式报告里其实很危险因为“p0”意味着“绝对不可能发生”这在统计上不成立。正确写法是“p 0.001”或“p 2.3e-16”。更重要的是理解“统计显著”不等于“实际显著”。一个20年NDVI趋势Sen斜率只有每年0.0005p值却高达0.001统计上显著因为序列非常平滑、噪声极小。但这个量级在生态学上可能完全没有意义。我每次写结论都会强制要求自己把“斜率数值的单位”和“实际物理意义”同时写出来。比如写成“NDVI平均每年上升0.002约等于20年累计0.04占研究区多年平均NDVI的2.1%”这样读者才有一个直观尺度。否则一个苍白的小数点数字既无法支撑结论也无法横向对比。5. 我的实操体会与几条进阶建议这套流程我用了很多年最大的感受是SenMK不是“一键出图”的黑盒它的门槛在原理理解、数据清洗和情境判断。同一份数据年度序列和月度序列用错方法结论可能完全相反同一份序列是否做自相关修正p值可能从0.04跳到0.20。所以每次跑完结果我都建议你花十分钟做三件事看一眼Sen斜率是否在物理合理范围内看一眼没通过显著性的区域是不是样本量特别少再看一眼空间分布是否出现了“地图上突兀的点状异常”——如果某个像素和周围大片区域趋势明显不一致先别急着下结论回到原始数据去查那个点是不是掩膜没做干净。如果数据里明显存在突变比如某一年之后趋势方向改变了我会在SenMK之外再叠加一个Pettitt突变检验或者分段SenMK把“整体趋势”拆成“前半段”和“后半段”分别描述。这能很好地弥补MK只检测单调趋势的短板尤其在分析政策干预、极端气候事件影响时特别有用。还有一点如果你处理的是长时间序列且要做空间对比建议统一设置随机种子以保证可复现如果再配上有效的样本量栅格和p值栅格你的论文审稿人或项目验收方基本找不到“统计方法交代不清楚”的指责点。最后直接给你一条经验趋势分析只是第一步结论要落在“为什么变”上。SenMK告诉你哪里变了、变了多少、变化可不可信但它不告诉你原因。想做归因还需要结合土地利用、气候因子、人类活动等空间数据做进一步分析那才是另一个有意思的开始。
Sen斜率与Mann-Kendall检验：时间序列稳健趋势分析实战指南
搞地学、遥感、水文、气象数据分析的老哥老姐们SenMK趋势分析这组词估计十有八九都熟。它基本是“时间序列趋势检测”里的默认组合了Sen负责估算变化速率MK负责判断趋势显著性两个一配合既能告诉你“涨了还是跌了”又能告诉你“这个涨跌能不能信”。我自己日常处理的不是NDVI就是降水、土壤湿度凡是论文里写“分析近二十年演变趋势”的基本都用这套方法交差。这篇文章就把我从原理理解、数据准备到Python跑批的完整流程连同踩过的坑一起拆开给刚开始接触趋势分析的学生也给每一个被长序列数据折磨的工程师。先给个结论这套方法不挑数据分布、对离群值不敏感、实现成本低属于那种“你不需要懂太多数学也能用对”的工具。但它也有边界比如对自相关数据、季节数据、短序列数据直接套用会出大问题。下面我会先把原理讲透再给完整代码和参数选择思路最后把高频问题列成排查清单照着做基本能避掉九成以上的坑。1. 为什么是“SenMK”这对组合解决什么问题1.1 常规线性回归的尴尬被离群值牵着走拿到一条时间序列多数人第一反应是做一元线性回归看斜率是否显著。这在教科书上没问题真实数据却常常不配合某一年降水特别大、某一年NDVI因为传感器问题突然跳变、或者源数据本身就带着几个极端噪声点。最小二乘回归对这类离群值非常敏感一个异常点就足以把斜率从正拽成负或者从显著拉成不显著。我印象很深的一次教训某站点20年降水序列我做线性回归得到的趋势是“每年增加约3.1毫米p0.02”看上去很漂亮。后来回查原始记录发现第7年有一次极端特大暴雨数值比正常年份高出好几倍。把这个点去掉再算回归斜率直接变成负的。一个点翻转整个结论这显然不合理。后来又用Sen斜率一算数字稳定在“每年下降约0.8毫米”因为Sen拿的是所有点对斜率的中位数那个极端暴雨点最多影响它参与的少数几个点对不足以撼动中位数。这也是为什么在真实项目中我宁可把线性回归当作初筛最终结论一律以SenMK为准。如果你不理解中位数的意义可以这样类比一群人里混进一个“年薪百亿”的算平均收入瞬间被拉爆但算收入中位数大家的生活水平几乎不受影响。Sen斜率就是“收入中位数版”的趋势估计它不追求完美拟合只求给一个不被异常值绑架的稳健方向。1.2 Sen斜率估计把所有点对斜率排个队Sen斜率也叫Theil-Sen估计计算起来非常朴素对于一条长度为n的序列任意取两个时间点i和ji 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分所有分数加总得到统计量SS Σ 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)(2n5) - Σ_t t_i(t_i-1)(2t_i5)] / 18后一项是“并列值修正”。如果数据里有大量完全相等的点比如降水经常出现0.0毫米的月份这个修正必须做否则方差被低估p值虚小容易得出“假显著”结论。Z值近似服从标准正态分布所以双侧检验p值就是2倍的尾部概率p 2×(1 - Φ(|Z|))。p小于0.05我们就说趋势显著小于0.01说极显著。MK检验不要求数据正态不怕离群值对大多数单调趋势的检测都很有效所以成了环境变化研究里的常青树。1.4 为什么两者总是成对出现Sen和MK能经常合体是因为它们优势互补、方法论基因一致。工具回答的问题输出Sen斜率变化速率是多快方向是增还是减斜率β单位是“每年多少量”MK检验这个趋势是否统计显著Z值、p值、Kendall TauSen算出“每年增加0.003的NDVI”MK负责告诉你“p0.42这个增加很可能是噪声”。反过来MK说“z2.31p0.02存在显著趋势”Sen则告诉你是每年1毫米还是每年50毫米。只用一个都容易误读一个“显著”但没有量级一个“有量级”但不可信。在实际项目中我还会顺带输出MK计算过程中得到的Kendall Tau它介于-1和1之间表示趋势的强度。写报告时表格里放“Sen斜率、p值、Tau值”三列基本就齐活了。2. 动手前的准备数据形态、工具选型与几项关键约定2.1 数据应该是什么样一维时间序列和多波段栅格SenMK的输入按项目场景分两种。第一种是单点时间序列比如某个气象站的逐年降水量或者某个采样点的水质指标一维数组长度为年份数。第二种是空间栅格比如MODIS NDVI的年合成产品20年就是20个波段每个波段代表一年的空间分布合成一个多波段GeoTIFF。处理栅格时本质上是把每个像素对应的“时间维数组”单独拎出来做一次SenMK再把结果写回成一个单波段栅格。无论哪种形态有一个前提必须搞清楚时间分辨率是什么。年度数据用普通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 matplotlibpymannkendall这个库非常完整原始MK、季节MK、自相关修正MK、预白化MK全都有函数名也直观。如果你不用PythonR语言里也有zyp和trend包功能类似但代码风格不同。这里我只讲Python路径毕竟栅格批处理这一块Python生态确实更顺手。2.3 两个常用实现路径手搓公式与现成函数做趋势分析最忌讳的就是黑盒使用。我的习惯是先用现成函数快速出结果再手写一遍核心公式验证理解。项目紧的时候直接调库项目松的时候至少把S统计量和Sen斜率自己算一遍这样出了意外结果才知道往哪个方向排查。pymannkendall最核心的接口是original_test它返回一个命名元组里面包含trendincreasing/decreasing/no trend、h是否拒绝原假设、pp值、zZ统计量、TauKendall Tau、sS统计量、var_sS的方差、slopeSen斜率、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_valid15): series np.asarray(series, dtypefloat) 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, dtypefloat) 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)你可能会担心这个双重循环太慢。确实n20时有190个点对n50时则上升到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(fSen斜率: {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, dtypefloat) 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(fS {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 Tauoriginal_test返回的slope就是第3.2节手写函数的结果两者应该对得上。注意result.trend对“是否拒绝原假设”的判断使用的是0.05显著水平如果你想要0.1或者0.01别直接用trend字段改用p值自己做判断。3.4 第四步批量处理栅格逐像素确实慢但别怕真正让SenMK项目从“练习”变成“实战”的是空间栅格批处理。假设你有一个NDVI年合成多波段GeoTIFFshape是(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(count1, dtypefloat32, nodatanp.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(fprocessed {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值最终画图时还要再分类。我常用的分级规则如下分级判定条件含义极显著上升Sen 0 且 p 0.01趋势非常可靠显著上升Sen 0 且 0.01 ≤ p 0.05趋势可靠弱显著上升Sen 0 且 0.05 ≤ p 0.1仅作参考趋势不显著p ≥ 0.1不能判断有趋势弱显著下降Sen 0 且 0.05 ≤ p 0.1仅作参考显著下降Sen 0 且 0.01 ≤ p 0.05趋势可靠极显著下降Sen 0 且 p 0.01趋势非常可靠出图时我习惯用暖色表示上升、冷色表示下降不显著区域统一用浅灰或透明。如果你处理的指标有物理范围比如NDVI理论范围是-1到1别忘了在色带数值上标注清楚单位比如“Sen斜率: 0.001/年”避免读者误读成“每年增加1%”。另外我还会生成一张“有效样本量”栅格当作辅助诊断图。当某些区域趋势看上去很奇怪时第一件事就是查那里的有效年份数是不是太少。很多时候一个显著性“突变区”其实是缺测造成的假象有了样本量栅格这种问题一目了然。4. 问题排查与避坑实录4.1 有效样本不足和NaN趋势结果不可信这是栅格处理中最常见的坑。MODIS NDVI常因为云覆盖导致某些年份某些像素为无效值降水数据遇到传感器故障也可能缺测。如果缺失年份太多MK检验的方差估计会失真。举例来说20年数据里只有8年有效序列长度和“自由度”都不足即使真的存在趋势也很难检测出来。我的做法是给每个像素设一个阈值比如有效样本≥15才返回结果低于阈值直接输出NaN并mask掉。这个阈值不是拍脑袋而是从检验功效角度考虑的n15时MK对中等强度的趋势尚能检测n8时基本什么都检不出来。如果你的时间序列本身就是短序列比如只有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里带自相关修正的方法。这个库提供了几个现成接口方法适用场景original_test独立同分布数据常规年度序列hamed_rao_modification_test方差修正版适合存在自相关的数据yue_wang_modification_test另一种方差修正方式trend_free_prewhitening_modification_test先去趋势再预白化适合正自相关明显的数据我的经验阈值是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)同理季度数据用period4周数据用period52。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(count1, dtypefloat32, nodatanp.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(windowwindow) # shape: (n_time, block_h, block_w) # 对block做像素级SenMK再把结果写回对应窗口计算过慢则取决于你点对循环的效率。如果年份数超过40每个像素点对数量至少780对全国范围内就是几十亿次比较。我的建议是核心函数用numba的njit装饰或者用numpy的向量化矩阵构建点对差分以牺牲部分内存换速度再或者用multiprocessing多进程并行。对单个项目而言先用几万像素小范围测试正确性再上全图避免程序跑了两小时才发现结果全是NaN。4.5 p值显示和解读别掉进“科学计数法”的坑MK输出的p值有时候是1e-16这种极小数字很多人偷懒写成“p 0.000”。这个习惯在正式报告里其实很危险因为“p0”意味着“绝对不可能发生”这在统计上不成立。正确写法是“p 0.001”或“p 2.3e-16”。更重要的是理解“统计显著”不等于“实际显著”。一个20年NDVI趋势Sen斜率只有每年0.0005p值却高达0.001统计上显著因为序列非常平滑、噪声极小。但这个量级在生态学上可能完全没有意义。我每次写结论都会强制要求自己把“斜率数值的单位”和“实际物理意义”同时写出来。比如写成“NDVI平均每年上升0.002约等于20年累计0.04占研究区多年平均NDVI的2.1%”这样读者才有一个直观尺度。否则一个苍白的小数点数字既无法支撑结论也无法横向对比。5. 我的实操体会与几条进阶建议这套流程我用了很多年最大的感受是SenMK不是“一键出图”的黑盒它的门槛在原理理解、数据清洗和情境判断。同一份数据年度序列和月度序列用错方法结论可能完全相反同一份序列是否做自相关修正p值可能从0.04跳到0.20。所以每次跑完结果我都建议你花十分钟做三件事看一眼Sen斜率是否在物理合理范围内看一眼没通过显著性的区域是不是样本量特别少再看一眼空间分布是否出现了“地图上突兀的点状异常”——如果某个像素和周围大片区域趋势明显不一致先别急着下结论回到原始数据去查那个点是不是掩膜没做干净。如果数据里明显存在突变比如某一年之后趋势方向改变了我会在SenMK之外再叠加一个Pettitt突变检验或者分段SenMK把“整体趋势”拆成“前半段”和“后半段”分别描述。这能很好地弥补MK只检测单调趋势的短板尤其在分析政策干预、极端气候事件影响时特别有用。还有一点如果你处理的是长时间序列且要做空间对比建议统一设置随机种子以保证可复现如果再配上有效的样本量栅格和p值栅格你的论文审稿人或项目验收方基本找不到“统计方法交代不清楚”的指责点。最后直接给你一条经验趋势分析只是第一步结论要落在“为什么变”上。SenMK告诉你哪里变了、变了多少、变化可不可信但它不告诉你原因。想做归因还需要结合土地利用、气候因子、人类活动等空间数据做进一步分析那才是另一个有意思的开始。