企业官网建设流程全解析

简介：面向气象、水文与环境数据分析学习者的Python实战脚本资源，聚焦季节尺度Mann-Kendall突变检测，用于识别具有周期性变化的时间序列中的趋势转折与显著突变点，适合课程实验、科研入门及方法复现。压缩包共2个文件，包含1个Python脚本和1个Excel数据文件，整体约11KB，脚本与数据配套，便于按步骤阅读、运行和核查结果。已有200人学习，说明该主题在相关分析场景中具有一定关注度。读者可从中获得M-K非参数检验的实现思路、季节性成分与趋势成分的分解方法、SPEI干旱指数数据的处理示例，以及突变点定位和结果可视化的参考写法，便于快速理解季节尺度突变检测的代码组织方式，并迁移到自己的时间序列分析任务中。

1. 从一个“整体显著”的假象说起

年尺度的 M-K 检验跑完，Z 值 -2.31，p=0.021，看上去是个干净利落的显著下降。可同一份月数据按季节拆开重新跑，夏季 Z 只有 -1.1，秋季掉到 -2.8，年尺度上那个“整体显著”其实全靠秋季一个季节撑着，春季甚至在往反方向走。这正是基于 Python 实现的季节尺度 M-K 突变检测要处理的事：先把月序列按季节聚合，再对每个季节单独做趋势与突变检验，最后合成一个能在季节之间横向比较的结论。它适合手里已经有月尺度长序列的人，降水、气温、径流、蒸散、NDVI、用电负荷都行，只要你不满足于“年总值一个 Z 值”的粗糙判断。数据量不大，一台普通笔记本几分钟就能跑完，改造成本也低。

2. 项目骨架与数据契约：一个 .py 脚本里其实藏着五层

单文件脚本最容易被写成“从上到下一条直线”，读三个月后自己都看不懂。这份实现按职责切成了五段：读数据、季节映射、M-K 统计量、突变扫描、出图落盘。每一段对外只暴露一两个函数，彼此之间靠 DataFrame 和 ndarray 传参，不共享全局变量。这样切的好处很实际——你可以只 import 统计量那一层去跑别的项目，也可以在季节映射那一层换成旬尺度、候尺度而完全不动后面的算法。

2.1 目录切分与模块职责

工程目录大致长这样，脚本名与主入口同名，方便直接python mk_seasonal.py：

mk_seasonal/ ├── mk_seasonal.py # 主入口：解析参数、串流程、写结果 ├── io_utils.py # 读入与清洗：长表 / 宽表两条路径 ├── season.py # 月 -> 季节映射、聚合、完整性校验 ├── mk_core.py # S 统计量、方差修正、Z 值、Sen 斜率 ├── mk_change.py # 顺序 M-K（UF/UB）与突变点判定 ├── report.py # 出图与 CSV/JSON 落盘 └── data/ └── monthly.csv # year,month,value 三列

mk_core.py是纯函数模块，不 import pandas，只吃 numpy 数组，好处是单元测试写起来极轻。season.py反过来只处理表格，不碰任何统计。主入口负责编排，把路径、显著性水平、季节方案这些参数从命令行接进来再往下传。我一般会把“会不会被复用”当作切分依据：统计算法一定会被复用，读数据的部分大概率不会，所以前者必须解耦。

2.2 输入数据契约：三列长表为主，宽表兜底

数据契约先定死：year、month、value三列，一行一个观测。这个约定看着朴素，但它挡住了后面 80% 的诡异 bug——重复月份、月份缺失、年份被读成字符串、value列里混进--之类的占位符。

# io_utils.py import pandas as pd REQUIRED = ("year", "month", "value") def load_monthly(path, value_col="value", na_policy="drop"): """读入月尺度序列，统一成 year/month/value 三列长表""" df = pd.read_csv(path) if value_col != "value": df = df.rename(columns={value_col: "value"}) missing = [c for c in REQUIRED if c not in df.columns] if missing: raise ValueError(f"缺少必需列: {missing}") df["year"] = df["year"].astype(int) df["month"] = df["month"].astype(int) # errors="coerce" 把 "--"、"NA"、空串统一变成 NaN，避免 object 列混进算法 df["value"] = pd.to_numeric(df["value"], errors="coerce") if na_policy == "drop": df = df.dropna(subset=["value"]) # 同一 (year, month) 出现多次时保留最后一条，防止合并表带来的重复 df = df.drop_duplicates(subset=["year", "month"], keep="last") return df.sort_values(["year", "month"]).reset_index(drop=True)

na_policy特意留了开关：drop用于趋势检验，缺测月份直接剔除；如果后面要画 UF/UB 曲线，缺测必须改成“打断序列”，不能简单删掉，否则时间轴被压缩，突变点会整体偏移，这条在第 5 章展开。drop_duplicates的keep="last"是按“后写入的更新”这个常见约定来的，如果你的数据源是先写入为准，改成first。

很多气象站点导出的其实是宽表，一行一年、列名是 1 到 12。转换走melt，不用手写循环：

def load_wide(path): """兼容 year x 12 月的宽表：一行一年，列名 1..12 或 '1'..'12'""" df = pd.read_csv(path) months = [c for c in df.columns if str(c).strip() in {str(m) for m in range(1, 13)}] long = df.melt(id_vars=["year"], value_vars=months, var_name="month", value_name="value") long["month"] = long["month"].astype(int) return (long.dropna(subset=["value"]) .sort_values(["year", "month"]) .reset_index(drop=True))

这里有个细节：months用集合匹配而不是df.columns[1:]全取，是为了防止宽表里还混着annual、station_id这种额外列被误 melt 进去。宽表路径我一般只在数据源确实是宽表时才用，主流程一律走长表，因为长表对缺测的表达是显式的——少一行就是缺，而不是一个空单元格。

2.3 季节映射与运行入口

季节映射是整个项目里最容易翻车的一环，不是算法难，而是“冬季跨年”这件事。12 月属于当年的冬季还是次年的冬季，不同学科的习惯不一样。代码里把选择权交出来：

# season.py SEASON_NH = {3: "MAM", 4: "MAM", 5: "MAM", 6: "JJA", 7: "JJA", 8: "JJA", 9: "SON", 10: "SON", 11: "SON", 12: "DJF", 1: "DJF", 2: "DJF"} def to_season(df, scheme="quarter", winter_label="end"): """month -> 季节标签；winter_label='end' 表示冬季按结束年份命名""" df = df.copy() if scheme == "quarter": df["season"] = df["month"].map(SEASON_NH) elif scheme == "month": df["season"] = df["month"].map(lambda m: f"M{m:02d}") else: raise ValueError(f"未知季节方案: {scheme}") # 12 月归入次年的冬季，保证 DJF 三个月落在同一个 season_year 上 shift = ((df["month"] == 12) & (winter_label == "end")).astype(int) df["season_year"] = df["year"] + shift return df

winter_label="end"是最常见的做法：某年冬季 = 上年 12 月 + 当年 1、2 月，按结束年份记。如果改成"start"，则冬季按起始年份记。两种都能用，但全流程必须统一，包括后面画图的横轴。

聚合方式取决于变量类型，降水、径流用sum，气温、NDVI、水位用mean，这个不能靠猜：

def aggregate(df, how="sum"): agg = (df.groupby(["season_year", "season"])["value"] .agg(how).reset_index() .rename(columns={"season_year": "year"})) return agg def season_completeness(df, expect=3): """检查每个季节是否凑齐 expect 个月，返回不完整的那些""" cnt = df.groupby(["season_year", "season"]).size() return cnt[cnt != expect]

主入口用 argparse 暴露参数，几个关键项是必须给出口径的：

python mk_seasonal.py \ --input data/monthly.csv \ --col value \ --agg sum \ --seasons quarter \ --alpha 0.05 \ --out results/

--agg决定聚合函数，--seasons决定走四季还是 12 个月，--alpha同时控制 Z 的显著性判读和 UF/UB 图上的置信线。把这三项提到命令行，是因为统计分析里，“换个口径再跑一遍”是常态，写死在代码里每次都得改文件。

提示：跑之前先season_completeness()过一遍，返回空 Series 再往下走，能省掉后面一半的排查时间。

3. 核心算法拆解：S 统计量、方差修正与 Sen 斜率

M-K 检验本身不复杂，难的是一堆看似不起眼的修正项。S 统计量是所有配对比较的符号和，方差在存在并列值（ties）时必须修正，否则 Z 值会系统性偏大，显著性虚高。这三块写成三个纯函数，谁都能单独验证。

3.1 S 统计量与结值修正

对长度为 n 的序列，S 定义为所有 i<j 的符号函数之和。直观含义是“后出现的值比先出现的值大的配对数 减去 小的配对数”，S 为正说明整体在上升。

# mk_core.py import numpy as np def mk_s(x): """返回 (S, Var(S))；Var 已做结值修正""" x = np.asarray(x, dtype=float) n = x.size if n < 3: return 0.0, 0.0 s = 0.0 for i in range(n - 1): # 向量化：一次比较第 i 个与它之后的全部元素 s += np.sum(np.sign(x[i + 1:] - x[i])) # 结值修正：相同数值分成若干组，t_p 为第 p 组的个数 _, counts = np.unique(x, return_counts=True) ties = np.sum(counts * (counts - 1) * (2 * counts + 5)) var_s = (n * (n - 1) * (2 * n + 5) - ties) / 18.0 return float(s), float(var_s)

np.unique返回的counts里，只出现一次的值 count 为 1，代入1*0*7=0，等于不贡献修正量，所以不需要额外过滤。这一步在降水量这种含大量 0 值的数据上影响非常明显——一个 40 年的月序列，零值可能占三成，不修正的话 Z 值能虚高 10% 以上。

注意：结值修正用的是“整个序列里相同值的分组”，不是“相邻相同值的分组”，别写成滚动比较。

3.2 方差、Z 值与显著性判读

Z 的计算有三段，S 大于 0 减 1，小于 0 加 1，等于 0 直接取 0。这个 ±1 是连续性修正，很多实现漏掉，样本量小于 20 时结果会差出可观的一截。

from scipy.stats import norm def mk_z(x, alpha=0.05): """返回 (Z, p, 是否显著)""" s, var_s = mk_s(x) if var_s <= 0: return 0.0, 1.0, False 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 - norm.cdf(abs(z))) return float(z), float(p), bool(p < alpha)

var_s <= 0这个分支不是防御性编程的摆设：当序列长度不足 3 或者全部数值相同时，方差会退化成 0，直接除会得到 inf 或 nan，后面画图就整个废掉。返回 Z=0、p=1 表示“无信息”，比抛异常更合适，因为批处理时你不想为一个站点中断整条流水线。

判读上要注意，Z 的绝对值大于 1.96 只说明“在 α=0.05 下趋势显著”，不说明趋势强弱。40 年序列和 10 年序列的 Z 值不可直接比较，样本量差异会放大前者。我一般会在结果表里同时给出 n，看的人自己心里有数。

3.3 Sen 斜率：比回归斜率更抗异常值的趋势幅度

Z 值告诉你“有没有趋势”，Sen 斜率告诉你“每年变多少”。它取所有配对斜率的中位数，对离群点几乎免疫，这是它比最小二乘回归更适合环境序列的原因。

def sen_slope(x, t=None): """Sen's slope：所有配对斜率的中位数""" x = np.asarray(x, dtype=float) n = x.size if n < 2: return float("nan") # 时间轴默认等距；有缺测时把真实年份传进来 t = np.arange(n, dtype=float) if t is None else np.asarray(t, dtype=float) i, j = np.triu_indices(n, k=1) # 上三角下标，正好覆盖所有 i<j slopes = (x[j] - x[i]) / (t[j] - t[i]) # 分母是时间差，不是下标差 return float(np.median(slopes))

np.triu_indices这一行省掉了双重循环，n=480（40 年月度）时配对数是 11 万级别，纯 Python 循环要跑几秒，向量化之后毫秒级。分母用t[j]-t[i]而不是j-i是有意为之：一旦序列中间有缺测年份被删掉，等距假设就破了，用下标差算出来的斜率单位是“每观测”，不是“每年”，量纲直接错。

结果单位是“每年”，如果原数据是月总量，Sen 斜率乘 12 才是年变化量。这个换算我一般写在报告层，不写进算法层，算法只负责给出最原始的数字。

4. 突变检测链路：UF/UB 曲线与季节维度扫描

趋势检验回答“变没变”，突变检测回答“哪一年变的”。M-K 的突变检测走的是顺序统计量路线：对序列的每个前缀算一次标准化统计量，得到 UF 曲线；对倒序序列同样处理再取负，得到 UB 曲线。两条曲线的交点就是候选突变点。季节尺度上，则是把这条链路在每个季节上各跑一遍。

4.1 顺序 M-K：UF/UB 两条曲线怎么算出来

顺序 M-K 用的是未修正的经典方差k(k-1)(2k+5)/72，因为前缀长度 k 在变化，逐个前缀做结值修正代价大、收益小，标准做法就是不修正。

# mk_change.py import numpy as np def sequential_mk(x): """返回与 x 等长的 (UF, UB)，第 0 个位置为 0""" x = np.asarray(x, dtype=float) n = x.size def _uf(seq): m = seq.size out = np.zeros(m) for k in range(2, m + 1): s = 0.0 for i in range(k - 1): s += np.sum(np.sign(seq[i + 1:k] - seq[i])) e = k * (k - 1) / 4.0 # 期望 v = k * (k - 1) * (2 * k + 5) / 72.0 # 方差 out[k - 1] = (s - e) / np.sqrt(v) if v > 0 else 0.0 return out uf = _uf(x) ub = -_uf(x[::-1])[::-1] # 倒序算 UF，取负再翻回来对齐原时间轴 return uf, ub

ub那一行是最容易写错的地方，展开说：先对倒序序列算出一组 UF，这组值的时间方向是反的；取负号得到 UB 的符号语义；再[::-1]翻转回原时间顺序。少任何一个操作，画出来的 UB 都会跟 UF 完全对不上，交点位置也就全错。

交点判定不能只看“两条线交叉”，必须加两个约束：交点落在置信带内，且交叉之后 UF 持续超出置信线。

def find_change_points(years, uf, ub, crit=1.96, persist=3): """候选突变年：符号翻转 + 落在置信带内 + 之后持续越界""" d = uf - ub pts = [] for i in range(1, len(d)): if d[i - 1] * d[i] >= 0: continue if abs(uf[i]) >= crit or abs(ub[i]) >= crit: continue tail = uf[i:i + persist] if len(tail) == persist and np.all(np.abs(tail) > crit): pts.append(int(years[i])) return pts

persist=3是我常用的经验值，意思是交点之后至少 3 个时间单位都待在置信带外。不设这个约束的话，一条噪声大的序列能报出五六个“突变点”，做报告的没法解释。

4.2 季节维度扫描：逐季 UF/UB 与合成检验

季节尺度有两种玩法，用途不同。第一种是逐季跑 UF/UB，看每个季节各自在哪一年突变，出 2×2 或 3×4 的面板图，适合找“夏季提前、秋季推迟”这类物候信号。第二种是合成检验，把各季节的 Z 合并成一个总体结论，回答“整条序列到底有没有显著趋势”。

合成用 van Belle-Hughes 的卡方分解，好处是不要求各季节样本量相等：

from scipy.stats import chi2 def seasonal_mk(df, alpha=0.05): """df 需含 season / year / value 三列，返回逐季结果与合成统计量""" rows = [] for season, g in df.groupby("season"): g = g.sort_values("year") z, p, sig = mk_z(g["value"].to_numpy(), alpha) rows.append({ "season": season, "n": len(g), "Z": z, "p": p, "significant": sig, "sen_slope": sen_slope(g["value"].to_numpy(), g["year"].to_numpy()), }) Z = np.array([r["Z"] for r in rows]) m = Z.size chi_total = float(np.sum(Z ** 2)) # 总体趋势 zbar = float(np.mean(Z)) chi_hom = chi_total - m * zbar ** 2 # 季节间是否一致 summary = { "chi2_total": chi_total, "df_total": m, "p_total": float(1 - chi2.cdf(chi_total, m)), "chi2_homogeneity": chi_hom, "df_hom": m - 1, "p_hom": float(1 - chi2.cdf(chi_hom, m - 1)), } return rows, summary

chi2_total检验“所有季节合起来有没有趋势”，自由度是季节数 m。chi_hom检验“各季节的趋势是否一致”，自由度 m-1。后者才是季节尺度分析的价值所在——如果p_hom小于 0.05，说明季节之间方向或强度不一致，这时候报一个“总体显著下降”就是误导，必须逐季报。

sen_slope传了真实年份作为时间轴，所以每个季节的斜率单位都是“每年”，四季之间可以直接比大小。这一点比很多人直接把季节序号当 x 轴要靠谱得多。

4.3 出图与结果落盘

图要能直接进报告，所以置信线、零线、图例、中文标注都得齐。matplotlib 的中文字体在无桌面环境里经常掉字，落盘前显式指定字体族，或者干脆全部用英文标签，我倾向前者。

# report.py import matplotlib matplotlib.use("Agg") # 无桌面环境必须，否则可能卡住 import matplotlib.pyplot as plt from scipy.stats import norm plt.rcParams["font.sans-serif"] = ["SimHei", "DejaVu Sans"] plt.rcParams["axes.unicode_minus"] = False # 负号显示成方块就是这行没设 def plot_ufub(years, uf, ub, ax=None, title="", alpha=0.05): crit = norm.ppf(1 - alpha / 2) ax = ax or plt.gca() ax.plot(years, uf, color="#c0392b", lw=1.6, label="UF") ax.plot(years, ub, color="#2c6fbb", lw=1.6, label="UB") ax.axhline(crit, ls="--", lw=1, color="gray") ax.axhline(-crit, ls="--", lw=1, color="gray") ax.axhline(0, lw=0.8, color="black") ax.set_xlabel("year"); ax.set_ylabel("统计量") ax.set_title(title, fontsize=10) ax.legend(loc="upper left", fontsize=8) return ax

matplotlib.use("Agg")这行放在 import pyplot 之前，顺序错了不生效。服务器上跑批时忘了它，进程会挂在后端初始化上，看起来像死循环。

四季面板图用 subplots 拼，横轴统一范围，方便横向对比：

def plot_seasons(season_data): """season_data: {季节: (years, uf, ub)}""" fig, axes = plt.subplots(2, 2, figsize=(10, 7), sharex=True) for ax, (name, (y, uf, ub)) in zip(axes.ravel(), season_data.items()): plot_ufub(y, uf, ub, ax=ax, title=name) fig.tight_layout() fig.savefig("results/ufub_seasons.png", dpi=200) return fig

落盘除了图，还要有能进一步处理的表。逐季结果写成 CSV，合成统计量写成 JSON，突变年单独一列，方便下游直接读：

import json, pandas as pd def dump(rows, summary, out_dir="results"): pd.DataFrame(rows).to_csv(f"{out_dir}/seasonal_mk.csv", index=False, encoding="utf-8-sig") with open(f"{out_dir}/summary.json", "w", encoding="utf-8") as f: json.dump(summary, f, ensure_ascii=False, indent=2)

CSV 用utf-8-sig是为了让表格软件打开不乱码，这个细节不值得省。

5. 避坑与排查：五个让突变点悄悄漂移的细节

这套流程跑通不难，跑出可信的结果难。下面五条都是我在实际数据上撞过的，每条按现象、原因、解决写，遇到对得上的直接照方抓药。

5.1 冬季 Z 值和年值对不上

现象：年尺度 Z 是 -2.3，四季里冬季 Z 却是 +0.4，符号完全相反，怎么看都不合理。

原因：12 月被算进了当年的冬季，导致某个“冬季”实际只有 12 月一个月的数据，另一个冬季有 1、2 月两个月，季节序列被切碎了，长度和年份都对不上。

解决：确认winter_label口径统一，用season_completeness()检查每个季节是否凑齐 3 个月，不完整的季节要么剔除要么标记。冬季序列的第一年和最后一年通常是残缺的，直接丢掉这两年比强行补齐更诚实。

5.2 Z 值比预期大一圈，显著性虚高

现象：同样是 40 年月度降水，别人算出来 Z=-1.8，你算出来 -2.4，卡在显著性边界两侧。

原因：要么漏了结值修正，要么用了未修正的方差；如果数据里还有明显自相关（月降水常年有滞后相关），独立同分布假设本身也不成立。

解决：先确认mk_s里的 ties 修正生效，拿一组含大量零值的序列手算验一遍。自相关的问题靠预白化处理，先拟合 AR(1) 去掉滞后相关再检验，代价是趋势会被削弱一点，但结论更稳。

5.3 UF/UB 交叉了四五次，不知道报哪个

现象：图上两条线缠成一团，交点一堆，每个交点后面 UF 都没能持续越界。

原因：序列噪声大或者长度短，find_change_points只做了符号翻转判定，没加持续性约束。

解决：加上persist参数，要求交点后连续若干个时间单位都在置信带外；同时把crit从 1.96 提到 2.58（α=0.01）再跑一遍，两次都保留下来的交点才写进结论。只报一个主突变点，其余作为次要候选列在附录。

5.4 删掉缺测月份后，突变年整体前移

现象：某站点有几年缺测，清洗时dropna一删，UF/UB 图上的交点比邻站早了三年。

原因：顺序 M-K 假设时间轴等距，删掉一行之后后续所有点的时间坐标都被压缩，累积统计量的对齐基准变了。

解决：缺测不要删，改成“打断”——在缺测处把序列切成若干连续段，每段单独跑，或者把缺测位置填 NaN 并在画图时断线。趋势检验可以删，突变检测不能删，这两条要分开处理。

5.5 合成卡方显著，但逐季看全都边缘

现象：chi2_total的 p 值 0.01，看着很显著，可逐季 Z 全在 -1.5 到 -1.9 之间，没有一个是显著的。

原因：这是合成检验的正常表现，多个弱信号叠加确实能显著。问题出在解读上，把它说成“所有季节都显著下降”就错了。

解决：报告里必须同时给出p_total和p_hom，以及逐季 Z 和 p。趋势方向一致但都不显著时，表述成“整体呈下降趋势，各季节均未达显著水平”才是准确说法。

提示：每次跑完先看p_hom，再看逐季 Z。p_hom小于 0.05 就别再谈总体结论了，直接讲季节差异。

6. 二次开发：把它改成能跑批、能出报告的工具

单个站点跑通之后，下一步基本都是“我有一百个站点的文件”。最省事的改法不是加循环，而是把主流程抽成一个函数，main只负责遍历和汇总。

# batch.py from pathlib import Path from io_utils import load_monthly from season import to_season, aggregate, season_completeness from mk_change import sequential_mk, find_change_points from report import dump def run_one(csv_path, agg="sum", alpha=0.05): df = load_monthly(csv_path) df = to_season(df, scheme="quarter", winter_label="end") bad = season_completeness(df) df_season = aggregate(df, how=agg) rows, summary = seasonal_mk(df_season, alpha) change_years = {} for season, g in df_season.groupby("season"): g = g.sort_values("year") uf, ub = sequential_mk(g["value"].to_numpy()) change_years[season] = find_change_points(g["year"].to_numpy(), uf, ub) summary["change_years"] = change_years summary["incomplete_seasons"] = int(len(bad)) return rows, summary if __name__ == "__main__": all_rows = [] for p in Path("data/batch").glob("*.csv"): rows, summary = run_one(p) for r in rows: r["station"] = p.stem # 站点标识跟着结果一起走 all_rows.extend(rows) dump(all_rows, {}, out_dir="results")

station字段在循环里补，不要放进seasonal_mk内部，那会让统计函数依赖外部命名。返回incomplete_seasons的计数也是个实用设计，跑完一批先看哪些站点数据不完整，再决定要不要重跑。

几个常见的扩展方向和改法，列在下面方便对号入座：

验证改造是否正确的办法很朴素：拿一组人工构造的序列，前 20 年恒定、后 20 年线性下降，突变点必须落在第 21 年。这个自检用例我每次动算法层都会跑一遍，比对着真实数据反复调参快得多。

血泪经验是，季节尺度的坑几乎全在数据准备阶段，算法本身反而稳。从那以后我每次换新数据集，都强制先跑一遍完整性检查和那个 20+20 的人工序列自检，两关都过了才看真实结果。希望帮到你。

本文还有配套的精品资源，点击获取