做气候数据分析这些年,我打交道最多的区域之一就是华北平原。这片土地上的气温和降水序列,表面看起来是一根根平滑的曲线,真正拆开看,年代际波动的痕迹到处都是。前阵子我系统做了一次华北平原的年代际气候变化、趋势和未来情景,从数据清洗一直做到未来情景集合预估,中间踩的坑、悟出来的门道,值得好好整理一遍。这篇文章不打算写成学术报告体的说明书,而是按我实际操作的顺序,把为什么选这些方法、参数怎么定、哪些结果必须小心解释,一五一十讲清楚。无论你是刚接触气候统计的新手,还是已经被业务数据折磨过几轮的老朋友,相信都能从中找到能直接上手的部分。
那会儿接到任务的时候,我心里很清楚:这不是简单画几根折线图就能交差的事情。华北平原作为重要的粮食产区,气候变化的影响会直接传导到灌溉用水调配、冬小麦和夏玉米的种植安排,以及极端旱涝事件的应对。所以我在设计分析流程时,先给自己定了三条原则:数据必须经过严格质检和均一性检查;趋势和突变结论必须有显著性检验支撑;未来情景必须说清楚不确定性从哪里来,不能把模式输出直接当预言用。这三个原则贯穿了后面所有的步骤。
1. 项目整体设计与研究思路
1.1 为什么把目光锁定在“年代际”尺度
先说清楚一个特别容易被绕晕的概念。日常天气预报关心的是未来几天,气候预测关注的是一个季节或一年的平均状态,而“年代际”一般指十年到三十年这个跨度的变化。它在时间轴上恰好夹在年际波动和长期趋势之间,属于“背景场”级别的信息。
我习惯用一个比喻来解释这种关系:年际变化就像股票每天的上下跳动,噪声非常大;而年代际变化相当于看它的月K线和年K线,把短期波动磨平之后,真正的方向才会显现出来。华北平原的年降水量本身就有很强的年际变率,丰水年能比枯水年高出一倍还不止。如果不把时间尺度拉长,很容易被某一年或某几年的极端情况带偏判断。比如在连续几个旱年后,可能就有人得出“华北正在持续变干”的结论;可一旦放到三十年尺度上看,可能只是一段干期在年代际背景里的位置。
对水资源规划、农业结构布局来说,几十年尺度的“背景变化”比单年的极端值更有决策价值。水库调度、灌溉定额、种植制度调整,这些动作的回报周期都是以十年为单位计算的。如果我只看最近五年的数据来做结论,基本上是在拿运气赌方向。
1.2 研究范围与要回答的问题
这个项目的空间范围很明确:华北平原,也就是我国东部北方那一大片冲积平原,基本算是一马平川的地形,海拔不高,非常有利于冷暖空气的交汇和干湿阶段的转换。为了照顾数据质量和代表性,我的分析同时做了两种口径:一种是区域平均序列,把所有站点先各自序列化,再做区域算术平均;另一种是保留站点的空间分布,看趋势在不同纬度和地形位置上的差异。
时间范围我选的是1961年到2023年。这个起点的好处在于,大部分观测站已经从人工观测转入相对规范的记录体系,数据长度也够算气候平均值和趋势检验。在气象业务里,足够长的基准期能有效减小气候态偏移带来的误差。
研究问题被拆成了四个:第一,华北平原的气温和降水在过去六十年里到底有没有显著变化;第二,变化过程中有没有明显的年代际转折点或突变;第三,干湿组合的变化趋势是什么,对农业生产意味着什么;第四,在不同排放情景下,未来三十年与本世纪中叶的气温和降水会走向哪里。这四个问题对应着后面的主要章节,也正好构成了一条完整的气候变化评估链条。
1.3 数据来源与工具准备的取舍
数据是整条分析链的基础,这里我踩过不少坑,尤其不能上来就拿一个数据集直接开跑。我用到的数据可以分为三类。
第一类是站点观测数据。国内的气象数据公开平台能申请到长序列的日值数据,包括气温、降水、相对湿度这些基本要素。站点数据的好处是观测相对可靠,缺点也很明显:站点分布不均,平原腹地密集,边缘山区稀疏,如果直接用单站结果代表整个区域,空间代表性会出问题。
第二类是网格化观测产品。这类产品已经做了插值和地形订正,站点点位的影响被缓解了很多,适合做区域平均。但它也有自己的毛病,不同产品的插值算法不同,在降水这种空间变率大的要素上,产品之间的数值差异可能超过趋势本身的量级。所以我做任何结论之前,都会至少拿两种产品做交叉验证。
第三类是未来情景数据,这里用的是国际主流气候模式比较计划第六阶段的产品。这套数据包含多个全球气候模式在不同共享社会经济路径下对未来的模拟结果。我对情景的理解是,它们并不是对未来的预言,更像是几组“如果这样发展下去会怎样”的推演。低排放、中等排放、高排放三个情景分别体现了不同的减排力度,后续分析会体现在结果应用上。
工具方面,我全程用Python完成。pandas处理表格和重采样,scipy做统计检验和线性回归,matplotlib和cartopy负责绘图。有些专门的统计方法不需要重新造轮子,比如Mann-Kendall检验,Python里有现成实现,我直接调包后再辅以自写代码核验关键参数的逻辑,这样既高效又能保证自己理解每一项输出的含义。
2. 数据预处理与关键指数构建
2.1 先把站点数据的“脏东西”挑干净
很多人拿到数据的第一反应是直接算平均值,这个习惯最好改一改。我在处理华北站点数据时,发现的问题主要有三类:降水出现负值、气温出现一秒内跳变十度以上的怪数据、以及长时间缺测。
负降水基本是仪器或者人工录入错误,直接剔除。气温跳变要谨慎,不能一看异常就删,得核对原始记录和相邻站点的同时段数据。比如某站某日气温从10℃瞬间跳到25℃再跳回来,而方圆五十公里内的其他站点都没有这个信号,那基本可以判定为错误。缺测处理最麻烦,我采用的策略是:如果当月缺测天数少于或等于5天,就用该月其余日数求平均作为月值替代;如果缺测太多,就把这个月标为缺失值,不做插补。整条序列中插补产生的数据如果超过总长度的5%,我会在报告里单独标注,因为这些人工填补的值可能携带额外的低频信号。
另一个容易忽略的问题是站点迁移。有些站点观测位置发生过变动,虽然公开数据里做过部分均一化,但直接拼接出的序列仍可能出现物理上说不通的断点。我处理时会先对每条站点序列做一次相对均匀性检验,发现断点后要么用断点前后的多年差值做修正,要么干脆把断点前段截掉,宁可用短一点的连续序列,也不给趋势分析埋雷。
2.2 构建真正有业务含义的指数
区域平均气温、年降水量这种最基础的指数当然要算,但它们不足以回答农业问题。我额外构建了一组关键指数。
气温侧,我算了年平均气温、夏季平均气温、冬季平均气温,还有一个对作物特别重要的指标:日平均气温稳定通过10℃的积温。这个概念可能有人不熟悉,简单说就是每个日平均气温超过10℃的日子,把超过部分的温度累加起来,反映的是热量供给。小麦和玉米的生长速度、熟期选择都和它直接相关。华北平原本来就是光热资源比较丰富的区域,积温如果持续增加,对冬小麦来说并不全是好事,因为暖冬会导致冬性品种的春化需求不满足,返青和抗寒都会受影响。
降水侧,除了年降水量,我还算了降水量距平百分率、夏季降水占比,以及标准化降水指数SPI。SPI的计算分两步:先用月降水序列拟合Gamma分布,再把累计概率转化成标准正态分布的分位数。SPI为负表示偏旱,为正表示偏湿,经过标准化处理后,不同时间尺度的干湿状态可以直接比较。这里要提醒的是,SPI只针对降水异常,没有考虑气温增高带来的蒸发需求放大。在气候变暖背景下,如果只看SPI,可能会低估实际旱情。要评估农业干旱,还得耦合潜在蒸散量,这种“暖干化”的问题在华北这样的半湿润区尤其需要注意。
2.3 关键代码:重采样与滑动平均的实现
数据准备部分的代码其实不复杂,但细节决定成败。我习惯把站点数据统一整理成“时间索引的月值表”,再用pandas重采样。下面是气温处理的核心代码,降水思路完全一致,只是聚合函数从平均变成求和。
import pandas as pd import numpy as np # 读取月值数据,时间列已经解析为datetime类型 df = pd.read_csv('华北区域_月值.csv', parse_dates=['time']) df = df.set_index('time') # 月均气温做年均,降水量做年总 annual_tas = df.groupby(df.index.year).mean() annual_pr = df.groupby(df.index.year).sum() # 对区域平均序列做5年滑动平均,突出年代际信息 sm_tas = annual_tas.mean(axis=1).rolling(window=5, center=True, min_periods=3).mean() sm_pr = annual_pr.mean(axis=1).rolling(window=5, center=True, min_periods=3).mean()这里把滑动平均的窗口参数设为5年,并且用中心化处理。中心化的意思是,2000年的滑动平均值实际上是1998到2002年的平均,而不是1995到1999年的平均。这个细节很多人不注意,直接决定一个转折点的位置判断。中心化滑动平均能减少相位偏移,让年代际阶段的位置更贴近真实。min_periods=3允许序列首尾有少量缺失时不至于完全丢光,但解读首尾若干年的滑动值时必须谨慎,因为那里的样本量不足,抖动会被放大。
3. 趋势、突变与年代际信号提取
3.1 线性趋势估计:一次回归未必能讲完整个故事
最传统的做法是用最小二乘回归求年值序列的线性趋势,斜率代表每年变化多少毫米或多少摄氏度,再把斜率乘以10换算成“每十年”的变化量,气象上这个单位用得最多。
from scipy import stats years = annual_tas.index tas_mean = annual_tas.mean(axis=1).values slope, intercept, r_value, p_value, std_err = stats.linregress(years, tas_mean) trend_per_decade = slope * 10 print(f'气温趋势:{trend_per_decade:.2f} ℃/10年,p值={p_value:.4f}')p值小于0.05才算在95%置信度下显著。但我要强调一个实操里最常见的误解:p值显著只说明趋势在统计意义上不是零,不代表趋势一定很大或者一定有用;反过来,p值不显著也不代表序列没有变化。当降水年际方差很大的时候,即使存在真实的年代际变化,线性回归也未必能给出显著性,这个特征在华北降水分析里非常典型。所以我不会只用线性趋势下结论,而是搭配非参数检验和滑动平均共同判断。
3.2 Mann-Kendall检验与自相关陷阱
Mann-Kendall检验在气候水文领域几乎是标配,原因是它不要求数据服从正态分布,对缺失值也相对宽容,还能抵抗少量极值的影响。它的核心逻辑是把所有数据两两配对比较,统计上升对和下降对的差值,再通过与随机情况对比来判断整体趋势。
import pymannkendall as mk res_tas = mk.original_test(tas_mean) print(res_tas.trend, res_tas.p, res_tas.slope)拿到输出之后,我要泼一盆冷水:Mann-Kendall检验有一个容易翻车的核心假设,就是序列内部互相独立。但气象序列往往存在正自相关,也就是今年的状态和去年有关联,这种自相关会虚增检验的显著性,把本来没那么强的趋势判定成显著。解决的办法之一是预白化处理:先对原序列拟合一阶自回归模型,把趋势信息与自相关信息分离,再用处理后的序列做MK检验。
预白化和趋势项的处理顺序要小心。常见做法是先估计并剔除线性趋势,再对残差估计自回归系数,最后把“原始序列减去自回归部分”得到白化序列,再回到原序列做检验。这一步做错顺序,会把趋势本身也一起抹掉。我在初学阶段就犯过这样的错误,导致温度趋势被压制得看不见,后来反复对照相关文献才意识到是预白化顺序的问题。
从实际分析结果看,华北平原的气温序列MK检验几乎年年都是显著上升,这点没有太多争议;但降水序列的MK结果在不同站点和不同时期差异极大,很多站点不显著。这本身就是一种信息:降水的变化更多体现为振荡和阶段性位移,而不是单调趋势。
3.3 突变检验:找到转折点,而不是猜转折点
在年代际分析里,比趋势更重要的是找清楚气候状态的阶段转换。我用的方法是滑动t检验和Pettitt检验搭配使用。滑动t检验的思路是把序列从中间某个位置切开成前后两段,比较两段均值是否有显著差异,通过滑动不断移动切割点,找到差异最大的位置。Pettitt检验则是一种非参数的突变检验,不要求正态性,通过比较序列所有位置前后排名的差异来定位突变点。
滑动t检验有三个参数要事先定:子序列长度、滑动步长和置信水平。子序列长度不能太短,太短会被年际噪声干扰;也不能太长,太长会把若干个阶段揉在一起。我在华北的分析中试过5年、7年、10年几个窗口,最后结合降水阶段规律选了7年。这个选择没有绝对标准,本质上是分析者在“分辨率”和“稳定性”之间找平衡。建议不要只用一个参数跑一遍就下结论,把几个窗口的结果都画出来,看转折点是否稳定出现。如果换个窗口转折点就跑掉了,那说明突变证据太弱,只能叫阶段性波动。
用温度序列做滑动t检验时,通常会看到一个上世纪80年代末到90年代初的明显转折,统计检验能给出较高置信度。降水序列则会呈现出十年级别的干湿交替,比如一段时间偏旱、一段时间相对偏湿。这个现象本身比任何单一趋势都重要,因为它是理解华北水资源紧张的关键背景。
3.4 EEMD分解:从混合信号里剥离年代际分量
趋势检验告诉我们方向,但想要把年际信号、年代际信号和长期趋势分开看,还得用更精细的信号分解工具。我给项目补充了集合经验模态分解,也就是EEMD。这个方法把复杂序列分解成若干个固有模态函数加上一个残差项,分别对应从高频到低频的不同尺度信号。
EEMD的核心思想是给原始序列添加多组白噪声进行分解,然后对结果做集合平均,从而克服传统经验模态分解的模态混叠问题。它在气候分析里用得越来越多,但同样必须说一句:EEMD分解出的分量不一定都有物理意义,更多时候是“统计上有信号、物理上待解释”。我在项目里用它辅助确认气温序列里的年代际成分,但从不单独拿一个分量去做显著性论证,更不会把某一层分解结果直接当预测指标。
实际操作时,EEMD对序列两端非常敏感,边缘效应会让开头和末尾的分解结果明显不稳定。我的处理是先在两端做数据延拓,再用延拓序列分解,最后截取中间段作为有效结果。碰到这些细节,最好能把分解结果和原始序列叠在一张图上检查,眼见为实。
4. 结果怎么读:气温、降水与未来情景的关键判断
4.1 华北平原的主要变化特征
把全部方法跑完之后,我梳理出了几个比较稳的结果特征。气温方面,华北平原的整体增暖趋势比较明确。就我手头样本和公开历史数据交叉验证的结论来说,年均气温大致呈现每十年升高0.2到0.3℃的量级,冬季增温幅度略高于夏季,这一点直接体现在冬小麦越冬条件的变化上。气温变化不是线性均匀的,中间明显有跳跃式的上升阶段,大约在上世纪80年代末到90年代初完成了一次气候态抬升。
降水方面,年度降水量没有出现全区域一致的单调下降趋势,更多的是十年尺度的干湿摆动。我做出的一张区域平均降水序列图上,20世纪60年代相对偏多,中后期经历了一段时间的少雨阶段,90年代末到本世纪初期又整体偏干,最近十几年则出现了波动性回升。看整体趋势很像一条水平线穿过大幅震荡的信号,但拆开阶段看,水资源紧张程度的差别非常大。
| 要素 | 常用统计量 | 过去几十年常见现象 |
|---|---|---|
| 年平均气温 | 线性趋势 / MK检验 | 整体上升,冬季升温更明显,有阶段跃升 |
| 年降水量 | 线性趋势 / 滑动平均 | 线性趋势多不显著,年代际干湿交替显著 |
| ≥10℃积温 | 积温累计 | 积温增加,生长季热量条件改善 |
| 夏季降水占比 | 季节统计 | 年内集中度明显,极端强降水事件增多 |
我在做图时习惯把三样东西叠在同一张画布上:原始年值用浅色细线,5年滑动平均用粗线,线性趋势线用虚线。这样的组合图能同时回答三个问题:这些年到底波动多大、阶段背景大概往哪个方向走、线性趋势是否有参考意义。很多只发年值曲线的报告,一眼看过去全是锯齿,不容易形成判断。建议每张趋势图都至少把滑动平均那条粗线画出来。
4.2 干湿组合变化比单一降水更重要
真正让我感到需要警惕的,不是气温或者降水单独的变化,而是它们的组合。华北平原的气温在上升,潜在蒸散发随之增加;降水却没有呈现出同步的显著增补。这样一来,同样的降水量对应的实际干旱风险会抬高,大气的“需水”变强了,“供水”却没跟上。
这类现象在农业上有一个很直接的后果:灌溉需求量增加。冬小麦生长季正处于春季和初夏,如果这期间降水偏少,叠加温度升高带来的蒸散加大,土壤水分亏损就会更严重。即使年降水总量没有明显变化,季节分配的变化也会对作物需水关键期造成影响。我特意算了生长季内的降水集中度,结果发现华北降水在时间分布上变得更不均匀,汛期集中度在上升,非汛期降水贡献在下降。这意味着“下雨的时候用不完,不下的时候不够用”的错配风险在上升。
除去平均状态,极端事件更值得单独看。暖冬会让小麦提前返青,种植敏感性增加;一旦春季再出现冷空气,冻害风险反而上升。这些影响从来不是靠一条年均温度曲线能暴露出来的,而是藏在季节尺度的温差和降水分配里。SPI指数在连续若干年负值的时段里,会清晰显示极端干旱年组的聚集现象,这也是年代际信号最直观的体现。所以做干湿分析时,别只盯着年均值,把季节尺度的组合变化画出来,往往能看到比平均状态更危险的东西。
4.3 未来情景:千万别把模式输出当精确预测
未来情景分析是所有环节里最容易产生误解的地方。全球气候模式的空间分辨率通常比较低,直接拿来代表华北平原某几个站点,结果偏差会大到没有参考意义。我的处理流程是:先把各模式的格点数据双线性插值到统一网格,再裁剪出华北平原区域;计算每个模式基准期的气候平均值,与同一时期观测值之间的偏差;将这个偏差作为系统性误差,在全时段上做订正。相当于给每个模式做了一次“校准”,然后再做区域平均和集合统计。
关于共享社会经济路径,我这里给出直白的解释:低排放情景代表积极减排、绿色发展得到强化的世界;中等排放情景代表延续当前的发展节奏,减排幅度有限;高排放情景代表政策基本不做重大调整,继续依赖化石燃料。基于这三组情景,模式集合大致给出的故事线是:无论哪种情景,本世纪中叶以前气温都大概率继续上升,差距主要体现在本世纪后半段。降水方面的信度比气温低,很多模式给华北降水增多,但增幅和空间分布的离散度非常大。因此我的报告里宁可展示10%到90%的集合波动范围,也不给一个孤零零的中位数预测。
未来情景的价值不在“预言”某年某月的雨量,而在帮决策者理解不同发展路径下风险范围有多大。华北平原的水资源规划如果只盯着历史数据做最优估计,很可能一上来就顶着上限;如果参考高排放情景的下限,则会早知道最坏情形有多糟。这才是情景分析真正在业务上起作用的地方。
5. 实操中的坑:常见问题与排查速查表
这部分我按自己经验列成了一张速查表,每一条都是实际踩过的跟头。
| 现象 | 可能原因 | 处理建议 |
|---|---|---|
| 气温序列出现明显断崖或突跳 | 站点迁移或仪器更换,均一性不足 | 先做序列均匀性检验;断点前后多年差值修正或截断使用 |
| 降水序列存在大量缺测 | 数据补齐方式不当会制造低频信号 | 缺测少于5天时用月内均值,多于5天标为缺失,插补比例超过5%要提示 |
| MK检验结果与直觉相反 | 序列自相关导致检验失效 | 先预白化再检验,顺序必须“先剔除趋势、再算自相关系数、最后做检验” |
| 同一套数据滑动平均视觉趋势明显但MK不显著 | 年际方差太大,线性趋势淹没在噪声里 | 改用滑动平均看阶段变化,同时给出分阶段统计,别把不显著直接说成无变化 |
| 不同网格产品算出同一时段降水趋势相反 | 插值方案差异导致区域代表性不同 | 至少用两种产品交叉验证,差异大时溯源站点数量与数据版本 |
| EEMD分解首尾两端分量异常摆动 | 边缘效应 | 两端延拓后再分解,分析时仅使用避开边界的中段 |
| 未来情景中某模式模拟的历史气温整体偏高 | 模式系统性偏差 | 用基准期偏差订正,不能直接使用原始格点值 |
除了表里的问题,我还想多说一个容易被低估的细节:滑动平均窗口长度的选择。我试过3年、5年、7年、11年四种窗口,同一套降水序列呈现出了完全不同的“阶段感”。窗口太长会把真实的干湿转折抹平,窗口太短又跟年际噪声纠缠不清。我的个人原则是:先用高频分量判断转折点候选位置,再用不同窗口交叉验证,最后结合物理常识确定一个能讲清楚故事的长度。在华北的降水分析里,5到7年的窗口都还算合理,超过10年就得反复掂量。
另一个值得重视的问题是站点密度不均匀带来的空间插值陷阱。华北平原核心区站点密集,边缘山地站点稀疏,如果用普通反距离加权插值,平原边缘会被内部站点“污染”,趋势的空间分布看起来过于平滑。我后来改用考虑地形约束的插值方法,并且把结果按站点进行对比,才算真正理顺了空间格局。空间插值没有银弹,关键是不管用什么方法,都要能回答“为什么这幅图值得相信”。
新手最容易犯的一个错误是拿高排放情景的某一个模式结果当作“预测值”写进摘要。其实单模式的极端结果放在集合里经常被拉回中位数附近,尤其是温度变化相对稳定,降水离散度极大。建议每次在报告里同时附上集合均值和分位区间,并且明确标注“不确定性范围”,这样至少不会让读者误以为未来可以准确到具体毫米数。
6. 我对这套流程的复盘与准备继续深入的几件事
整套分析做下来,我最深刻的体会是:气候变化的“显著”和“有用”是两回事。统计显著只说明信号不容易被噪声解释,但真正对农业、水利有意义的,是信号背后的物理机制和应对空间。华北平原的增温趋势非常稳定,这给热量资源带来的改变是确定的;但水资源怎么变化、极端旱涝怎么重新分布,目前的不确定性依然很大。与其纠结给降水预测“定一个数”,不如把多种情景下的风险区间画清楚,让决策者在知道下限和上限的情况下去准备。
还有一个方法论层面的体会:不要迷信任何一种单一统计方法。线性趋势、MK检验、突变检验、EEMD,每种工具都有自己的假设窗口和盲区,只有把它们叠在一起交叉验证,结论才站得住。我在不同阶段用同一批方法得出来过互相矛盾的表象,最后都是靠回到原始序列、回到站点分布、回到物理机制来解开的结。气候统计是一个需要反复自省的工作,很多时候数据会说谎,说谎的源头不是数据本身,而是分析者没有认真理解方法的前提条件。
后面我打算把这件事往两个方向继续推进。一是把作物生长模型接进来,把气温和降水变化直接翻译成冬小麦和夏玉米的产量与需水变化,让分析结果离田间更近一步。二是引入更高分辨率的区域气候模拟数据,专门研究极端降水和高温热浪在华北平原不同区域的演变细节。年代际分析不应该停在纸面上的统计游戏,它最终要能回答:明天要种什么,地要怎么浇,风险要用什么方式扛下来。这条路我会继续往下走,也欢迎你也把数据、方法和结论摊在桌面上一起聊。