6.1 综合案例：气象数据分析管线

6.1 综合案例：气象数据分析管线

本节摘要：气象数据分析管线是把全书模块串成一条完整流水线的综合案例：带缺失与噪声的站点气温数据，先经插值补全，再做平滑滤波，随后用傅里叶变换找出季节周期，用线性回归检验长期趋势，最后借助 KDTree 找到邻近站点并计算相关性，输出一份结论明确的报告。整条管线按函数组织、分步输出中间结果，是真实工程组织方式的完整示范。

你能学到什么

阅读完本节，你应当能够：

1. 独立完成"插值补全 → 平滑滤波 → 周期分析 → 趋势检验 → 空间相关 → 报告输出"六步管线；
2. 说清每一步为什么选那个模块、那个参数，以及不选的代价；
3. 用函数划分和中间结果输出组织一段科学计算代码；
4. 把固定随机种子、逐站循环、结果汇总这些工程习惯用到自己的数据上。

一、问题与直觉

假设你接手了一份三年的日平均气温数据，来自三个站点，每个站点大约 8% 的观测因为设备故障变成了空值，测量噪声也不小。领导给你四个问题：季节周期到底是不是一年？这些年有没有变暖趋势，显著吗？邻近站点的温度变化是不是同步的？最后要一份能直接读的报告。

这份数据我们直接生成模拟版本，好处有三：一是可复现，随机种子固定，你跑出来的数字和我们一致；二是每个环节的正确性可验证——比如我们明知真实周期是 365 天，就能检查 FFT 有没有找对；三是真实数据只需替换读取那一小段，其余管线原样可用。这是科学计算里常用的"先模拟后实战"策略，比拿一份来路不明的真实数据直接开干稳得多。

二、核心原理：把问题拆成六步

整条管线是一个串行流程，每一步的输出都是下一步的输入：

图：气象分析管线总览

先看数据生成与预处理部分。我们生成三个纬度的站点：北京、上海、广州，基线温度依次升高，各自叠加季节项、缓慢升温项和噪声：

import numpy as np
from scipy import interpolate, signal, fft, stats, spatial

rng = np.random.default_rng(42)

def generate_series(days, base_temp, rng):
    """生成带季节、趋势与噪声的日平均气温序列。"""
    seasonal = 15.0 * np.sin(2 * np.pi * days / 365.0 - 1.2)
    trend = 0.002 * days
    noise = rng.normal(0, 1.5, size=days.size)
    return base_temp + seasonal + trend + noise

def add_missing(series, ratio, rng):
    """随机挖掉一部分观测，模拟设备故障。"""
    out = series.copy()
    idx = rng.choice(series.size, int(series.size * ratio), replace=False)
    out[idx] = np.nan
    return out

days = np.arange(365 * 3, dtype=float)
names = ['北京', '上海', '广州']
base_temps = [12.0, 17.0, 22.0]

stations = {}
for name, base in zip(names, base_temps):
    raw = generate_series(days, base, rng)
    stations[name] = add_missing(raw, 0.08, rng)
​

这里把"生成"和"挖缺失"拆成两个函数，是刻意的：真实项目里，读取数据、清洗、分析各归各的函数，才能单独测试、单独替换。default_rng(42) 固定了随机源，整套流程可以复现。

2.1 第一步：插值补全缺失值

def fill_missing(days, series):
    """用三次样条插值补全缺失值。"""
    mask = np.isfinite(series)
    interp = interpolate.interp1d(days[mask], series[mask], kind='cubic')
    return interp(days)
​

为什么用三次样条而不是线性？线性插值在缺口两端会留下折角，后面做平滑和周期分析时这些不光滑点会泄漏能量；三次样条保证导数连续，曲线更自然。代价是样条可能在缺口处轻微过冲——对 8% 的随机缺失来说，两种方法的差别不大，我们选三次样条是给后续环节留余量。

有两个细节值得注意。第一，interp1d 要求 x 严格递增，days 满足；第二，默认情况下在数据范围外求值会直接报错，这是一种防呆设计——补缺失是在范围内求值，没问题；如果你真想外推，必须显式传 fill_value='extrapolate'，并且想清楚外推结果靠不靠谱。

2.2 第二步：平滑滤波

def smooth_series(series, window=15, polyorder=3):
    """Savitzky-Golay 平滑，压低高频噪声。"""
    return signal.savgol_filter(series, window, polyorder)
​

signal.savgol_filter 是滑动窗口内做多项式拟合的平滑器，比起简单移动平均，它的优势是保形：不会把峰值削平、不会引入相位滞后。窗口必须取奇数，阶数必须小于窗口长度。我们取窗口 15 天、阶数 3——窗口远小于 365 天的季节周期，所以不会伤到季节信号，只压掉逐日的高频抖动。这就是"参数要配着信号的时间尺度选"的活例子：滤波器参数不是随便填的，它隐含了你对信号尺度的判断。

2.3 第三步：FFT 周期分析

def find_main_period(days, series):
    """用 FFT 找主周期，返回周期天数。"""
    centered = series - np.mean(series)
    spectrum = fft.fft(centered)
    freqs = fft.fftfreq(days.size, d=1.0)
    positive = freqs > 0
    peak = np.argmax(np.abs(spectrum[positive]))
    return 1.0 / freqs[positive][peak]
​

这里有个必踩的坑：不先去均值，0 频分量会占据绝对峰值，把周期信号全盖住。fft 返回的是复数谱，我们取模找最大峰值；fftfreq 的第二个参数是采样间隔，按天记就是 1.0，这样算出的频率单位是"每周期多少次"，取倒数就是周期天数。我们的模拟数据长度是 1095 天，正好是 365 的整数倍，频谱泄漏很小；真实数据往往凑不齐整周期，峰值会变宽、出现旁瓣，届时要考虑加窗函数，这点在第 4 章讲过。

2.4 第四步：趋势拟合与显著性检验

def trend_test(days, series):
    """线性趋势拟合与显著性检验。"""
    return stats.linregress(days, series)
​

stats.linregress 一次返回斜率、截距、相关系数、p 值和斜率标准误。注意斜率单位是"度每天"，写报告时要换算成"度每十年"，乘以 3650 即可。p 值小于 0.05 且斜率为正，才能说"显著变暖"。气象上更严格的做法会考虑序列的自相关对 p 值的影响，这里我们用教科书式的最小二乘检验，够用且好解释。

2.5 第五步：空间相关性

coords = np.array([[39.9, 116.4], [31.2, 121.5], [23.1, 113.3]])
tree = spatial.cKDTree(coords)
dist, idx = tree.query(coords, k=2)   # 第二近的才是真正的最近邻
​

spatial.cKDTree 建一次树，查询就是对数复杂度。query(k=2) 是因为每个站点最近的"邻居"是自己，第二近的才是要找的邻近站点。真实项目里可能有几百个站点，这种数据结构才扛得住。

然后是相关性计算。这里藏着一个容易犯的错：直接用原始温度序列算相关，所有站点的相关系数都会接近 1，因为季节成分是同步的，谁和谁都像。正确的做法是去掉季节项，用异常序列（距平）算相关，这才反映"今天北京偏暖时，上海是不是也偏暖"这种同步性：

for i, name in enumerate(names):
    neighbor = names[idx[i, 1]]
    a = results[name]['smooth'] - results[name]['smooth'].mean()
    b = results[neighbor]['smooth'] - results[neighbor]['smooth'].mean()
    r, p = stats.pearsonr(a, b)
    print(f"{name} 与最近邻 {neighbor} 的异常序列相关: r={r:.3f}, p={p:.3g}")
​

2.6 主流程与报告输出

results = {}
for name in names:
    filled = fill_missing(days, stations[name])
    smooth = smooth_series(filled)
    period = find_main_period(days, smooth)
    trend = trend_test(days, smooth)
    results[name] = dict(smooth=smooth, period=period, trend=trend)
    print(f"{name}: 主周期 {period:.0f} 天，趋势 {trend.slope * 3650:.3f} 度每十年，p 值 {trend.pvalue:.3g}")

print("--- 分析报告 ---")
for name in names:
    tr = results[name]['trend']
    if tr.pvalue < 0.05 and tr.slope > 0:
        verdict = '显著变暖'
    else:
        verdict = '无明显趋势'
    print(f"{name}: 季节周期约 {results[name]['period']:.0f} 天；{verdict}，升温 {tr.slope * 3650:.3f} 度每十年")
​

把上面的函数和主流程按顺序拼在一起，就是一份完整可跑的管线程序。跑出来的结果应当是：三个站点主周期都接近 365 天，升温速率在 0.7 度每十年上下，p 值极小，邻近站点的异常序列相关在 0.9 以上——因为我们的模拟数据是同一个季节项加独立噪声，这个数字合理。换成真实数据，解读方式完全一样，数字会告诉你真实世界的答案。

三、工程实践要点

3.1 环节、模块与参数对照

这张表也是你自己写管线时的检查单：每个环节都能说清"用什么、为什么、回答什么问题"，说明你真正理解了这段代码。

3.2 中间结果输出与函数划分

真实的科学计算程序，一半的功夫花在"让人看得懂中间发生了什么"。我们的做法是：每个函数只做一件事，输入输出都是数组，可以单独喂数据测试；主流程每步打印关键数字——补全率、主周期、斜率、p 值、相关系数——程序跑到哪一步、哪一步结果可疑，一眼就能看出来。这比写完一坨再回头调试省事得多。

还有一道看不见的工序：入口校验。真实数据不会像模拟数据这么守规矩——站点名可能缺失、坐标可能写反、缺测比例可能高达三成。建议在主流程开头做三件小事：打印每个站点的缺测比例，超过阈值就报警；检查坐标数组的长度与站点列表是否一致；对插值结果做一次有限性检查，确认补全后没有出现 NaN。这三行检查在模拟数据上看似多余，换成真实数据就是救命稻草。

⚠️ 常见坑：FFT 前忘了去均值，0 频峰值会盖住一切周期信号，找出的"主周期"永远是无穷大。另外 interp1d 默认禁止外推，想外推必须显式声明，而且先想清楚外推结果有没有物理意义。

💡 关键直觉：整条管线的设计原则是"数组进、数组出"——每步的输入输出都是一维数组，任何一步都可以单独调试、单独替换实现。你以后写自己的分析管线，照这个原则拆函数，永远不会拆出无法测试的意大利面。

3.3 一步一验：管线要边走边看

管线越长，越要防"最后结果错但不知道错在哪"。我们的经验是三步一验：第一步，每步打印中间结果，主周期接近 365、p 值数量级合理，才继续往下走；第二步，把中间结果和原始数据画在一张图上，肉眼确认平滑没有削掉季节峰；第三步，先在小样本上跑通全流程，再放开全量数据。这三步花的时间，远小于最后排错花的时间。

3.4 管线设计的两个取舍

第一个取舍：用 linregress 还是更复杂的季节分解？如果趋势和季节耦合严重，可以先用 FFT 去掉季节分量再回归，或者用 signal 里的滤波把季节带滤掉。我们的数据季节项固定，直接回归没问题；真实数据若有年际变化，建议先分解再检验。第二个取舍：站点相关性用异常序列还是原始序列？前面说过，异常序列才反映同步波动，这是"算对了"和"算得热闹"的区别——原始序列的相关几乎必然高，但那不是你想回答的问题。

本节速览

• 管线六步：插值补全、平滑滤波、FFT 周期分析、趋势检验、KDTree 邻近查找加相关分析、报告输出，环环相扣。
• 插值选三次样条：保证导数连续，利于后续平滑与频谱分析；默认不越界求值是防呆设计。
• savgol 保形平滑：窗口必须奇数、阶数小于窗口，窗口尺度要远小于目标信号周期。
• FFT 先去均值：否则 0 频峰值盖过周期峰；采样间隔参数 d 决定频率单位。
• linregress 一条龙：斜率、p 值一次拿全；斜率单位换算成每十年再写进报告。
• 异常序列算相关：原始序列的相关被季节主导，距平序列才反映真实的同步性。
• 工程组织：函数输入输出都是数组、每步打印中间结果、固定随机种子，是管线可交付的底线。

管线能跑之后，下一个问题自然冒出来：站点增加到几百个、序列拉长到几十年，还能跑得动吗？下一节我们讲性能优化——先测再优化，向量化、编译、稀疏化轮番上场。