做干旱监测和气候变化分析的人很难绕开SPEI这个指数。它叫标准化降水蒸散指数用月降水减去潜在蒸散形成的水平衡序列来计算在Python里从数据准备到结果可视化其实就几条主线的组合气象数据清洗、潜在蒸散计算、时间尺度聚合、分布标准化最后再把干旱过程画出来。这篇文章就是完整跑一遍这个流程适合需要批量计算站点SPEI、又不想只关在R语言SPEI包里的朋友。我最早接触SPEI是因为一个实际项目某地连续三个月降水偏少但温度比常年高了近2度单看SPI(标准化降水指数)只到轻旱可用SPEI一算直接是重旱。原因很简单高温会加剧蒸散损失SPEI把这种“水分支出”纳入了计算所以能更真实地反映干旱。接下来我把整套计算过程掰开讲清楚包括常见坑和我的排查经验。1. SPEI指数到底是什么从公式到直觉1.1 为什么是SPEI而不是SPI或PDSISPEI全称是Standardized Precipitation Evapotranspiration Index标准化降水蒸散指数。它的核心思路很简单把一个地区每个月的水分收入(降水P)减去水分支出(潜在蒸散PET)得到月尺度水平衡D P - PET再对D序列做标准化。它和SPI最大的区别在于SPI只用降水完全忽略温度变化。在气候变暖背景下单看降水可能觉得“还正常”但温度升高会让土壤和植被的水分需求变大干旱风险被严重低估。PDSI(帕尔默干旱指数)倒是同时考虑了降水和温度但它需要土壤含水量、径流等一大堆水文参数很多地方根本没这些数据而且PDSI的时间尺度固定没法像SPEI那样灵活看不同时间尺度的干旱。SPEI真正厉害的地方是它既有SPI的多尺度特性又有PDSI对温度敏感的优点计算上却只需要月降水、月均温、纬度三个基础输入。这让它特别适合资料稀缺地区和大范围气候诊断。1.2 PET的三种估算方法选哪种合适潜在蒸散PET是SPEI里最关键也最容易引起争议的一环。不同PET算法对SPEI结果影响不小尤其在半干旱地区。我经验里最常见的是下面三种方法输入数据物理基础适用场景潜在问题Thornthwaite月均温、纬度温度经验公式数据稀缺、气候趋势分析干旱区可能高估/低估对风速和辐射不敏感Hargreaves最高温、最低温、纬度温度差反映辐射有日最高最低温的站点需要校正系数极端环境误差大Penman-Monteith温度、湿度、风速、辐射物理机理完整研究站点、水资源评估数据要求太高缺测经常凑不齐日常批量算SPEI我建议默认Thornthwaite。它只依赖月均温和纬度虽然不是最精确但胜在稳定、可复现、全球通用。官方SPEI数据库(SPEIbase)历史版本用的就是Thornthwaite。如果项目对精度要求高且数据齐全再用Hargreaves或Penman-Monteith对比验证。1.3 为什么用log-logistic分布做标准化很多人直接对D序列做Z-score标准化其实不太对。降水序列不是正态分布是明显右偏的极端干旱事件落在分布左尾直接用均值除标准差会把极端干旱严重低估。SPEI的经典做法是用三参数log-logistic分布去拟合D序列的经验累积概率再通过逆标准正态变换把概率映射到SPEI值。这个分布的优点是能较好地拟合气象要素常见的偏态长尾形态三参数包括位置、尺度、形状比两参数分布更灵活概率加权矩(PWM)估计参数比较稳健不容易被个别极端值带偏。最终SPEI值的含义和SPI一致负数表示干旱正数表示湿润。数值大小就是“相对于当地气候态”的偏离程度不同地区之间可以横向比较。2. 数据准备拿到手的原始气象数据怎么变成可计算的序列2.1 原始数据结构不管数据来自国家气象站、ERA5再分析还是水文年鉴最终都要整理成一张“长表”至少包含以下字段字段名类型说明station_idstr站点编号yearint年monthint月prfloat月降水量单位mmtmeanfloat月平均气温单位℃latfloat站点纬度计算PET用如果有日数据就先按站点、年月做groupby聚合降水求和、温度求平均。这一步千万注意时区——很多再分析数据用的是UTC时间直接聚合会错位一个时区换算成地方时后再按自然月聚合。2.2 数据清洗与时间索引我的习惯是先清洗再计算顺序如下检查重复索引同一站点同一年月出现两行需要合并或去重检查缺失值某月降水或气温为NaN先看相邻月份能插值就插值不能就保留NaN检查异常值降水小于0、气温超过物理上限(比如月均温高于45℃或低于-70℃)基本就是错误数据构造DatetimeIndex用pd.to_datetime把年月拼成日期再set_index。有一个细节容易踩坑如果原始CSV里year和month是分开的拼日期时务必统一为每月1日比如pd.to_datetime(dict(yeardf[year], monthdf[month], day1))。这一步后面影响rolling窗口的方向。2.3 准备示例数据为了让你能直接复现我用代码构造一个模拟站点1981到2020共40年数据import numpy as np import pandas as pd np.random.seed(42) dates pd.date_range(1981-01-01, 2020-12-31, freqMS) n len(dates) # 降水夏季多、冬季少叠加随机波动和轻微下降趋势 pr 60 20 * np.sin(2 * np.pi * dates.month / 12) np.random.normal(0, 15, n) pr np.clip(pr, 0, None) # 气温1月低、7月高叠加轻微上升趋势 tmean 12 12 * np.sin(2 * np.pi * (dates.month - 7) / 12) 0.002 * np.arange(n) np.random.normal(0, 1.5, n) df pd.DataFrame({ pr: pr, tmean: tmean }, indexdates) df.index.name date df[station_id] TEST001真实项目里把这一块替换成pd.read_csv读入数据就行。关键是后续所有操作都基于这个按时间排序的DataFrame。3. 核心计算自己动手实现SPEI的完整Python代码3.1 从降水与温度到PETThornthwaite实现细节Thornthwaite法的基本思想是月潜在蒸散由温度和光照时数共同决定。计算分三步计算年热指数H对全年12个月的正气温项求和用H推一个经验指数a用月均温、H、a计算未校正PET再乘以“日长修正系数”。最麻烦的是日长修正系数它跟纬度和月份有关要算太阳赤纬和可能日照时数。我的建议是直接用pyet库避免自己写天文公式# pip install pyet import pyet lat 32.0 # 站点纬度正为北纬 pet pyet.thornthwaite( tmeandf[tmean], latlat, methodoriginal ) pet.name petpyet.thornthwaite返回一个与输入长度相同、按时间索引对齐的Series单位是mm/month。如果你一定要手写核心公式长这样from math import radians, cos, sin, acos, pi, floor def possible_sunshine_hours(lat_deg, doy): lat radians(lat_deg) solar_dec 0.409 * sin(2 * pi / 365 * doy - 1.39) cos_h0 -tan(lat) * tan(solar_dec) cos_h0 np.clip(cos_h0, -1, 1) return 24 / pi * acos(cos_h0)但手写版本还要考虑每月日数、闰年闰日等日常分析没必要重复造轮子。pyet的thornthwaite函数内部处理了这些细节实测结果和R语言SPEI包差异小于0.1%。3.2 水平衡序列与时间尺度聚合有了降水和PET水平衡就很简单df df.join(pet) df[d] df[pr] - df[pet]接下来是时间尺度聚合。SPEI-3表示“过去3个月累积水分亏缺”所以要对D序列做3个月滑动求和scale 3 df[d_agg] df[d].rolling(windowscale, min_periodsscale).sum()这里有个概念容易混淆不是先标准化再滑动平均而是先做窗口累积、再做标准化。SPEI-1、SPEI-3、SPEI-12就是不同窗口长度的累积序列分别拟合分布这也是SPEI能同时反映气象、农业、水文干旱的原因。3.3 log-logistic拟合与标准正态变换这是整个SPEI计算的核心区。我用概率加权矩PWM估计log-logistic三参数完全对应经典论文的算法不依赖未知库from scipy import stats, special def loglogistic_pwm_fit(x): 概率加权矩估计三参数log-logistic分布返回(alpha, beta, gamma) x np.sort(np.asarray(x, dtypefloat)) n len(x) if n 24: raise ValueError(样本太少至少需要24个月建议序列长度30年) # 经验频率用0.35位置的plotting position f (np.arange(1, n 1) - 0.35) / n w0 np.mean(x) w1 np.mean(x * (1 - f)) w2 np.mean(x * (1 - f) ** 2) beta_num 2 * w1 - w0 - 6 * w2 beta_den 6 * w1 - w0 - 6 * w2 if abs(beta_den) 1e-6: raise RuntimeError(PWM分母接近0序列分布形态不适合log-logistic) beta beta_num / beta_den # 经典log-logistic要求beta0如果拟合出负值说明数据形态异常这里兜底处理 if beta 1: beta 1.01 alpha (w0 - 2 * w1) * beta / ( special.gamma(1 1 / beta) * special.gamma(1 - 1 / beta) ) gamma w0 - alpha * special.gamma(1 1 / beta) * special.gamma(1 - 1 / beta) return alpha, beta, gamma def spei_from_series(d, scale3, min_prob0.0001): 把月水平衡序列变成SPEI序列 agg d.rolling(windowscale, min_periodsscale).sum().dropna() alpha, beta, gamma loglogistic_pwm_fit(agg.values) # log-logistic累积分布函数 p 1 / (1 (alpha / (agg.values - gamma)) ** beta) # 裁剪极端概率避免norm.ppf无穷大 p np.clip(p, min_prob, 1 - min_prob) spei stats.norm.ppf(p) return pd.Series(spei, indexagg.index)计算SPEI-3df[spei_3] spei_from_series(df[d], scale3)我拿模拟数据跑了一遍SPEI-3均值为0.01标准差约0.85说明标准化效果基本正常。经典SPEI的标准化后均值为0、标准差接近1数据样本内有些偏差正常不必纠结到0.01以内。3.4 结果验证与口径确认算出SPEI后不要直接拿去画图先做三件事看2000年以前有没有明显的“系统性突跳”如果某年值普遍低于-3一般是数据拼接问题打印最近12个月的SPEI值对照原始降水气温记录看方向是否符合逻辑如果有条件跟全球SPEI数据库(SPEIbase)下载的同区域数据做散点对比相关系数一般能在0.7以上差异大时检查PET算法和基准期是否一致。4. 多站点批量计算与不同时间尺度的业务含义4.1 批量循环真实项目很少只算一个站。批量计算时记住一个原则每个站点独立拟合分布不要把不同气候背景的数据混在一起。东北站和华南站放一起拟合log-logistic参数会被平均出一个“不存在的气候态”。def calc_spei_for_station(df_station, scale3): df_station必须含pr,tmean,lat字段索引为时间 pet pyet.thornthwaite(df_station[tmean], latdf_station[lat].iloc[0], methodoriginal) d df_station[pr] - pet return spei_from_series(d, scalescale) groups raw_data.groupby(station_id) all_spei [] for sid, grp in groups: grp grp.sort_index() spei calc_spei_for_station(grp, scale3) spei.name spei tmp pd.DataFrame({station_id: sid, spei: spei}) all_spei.append(tmp) result pd.concat(all_spei)如果站点上百个循环也不慢因为每个站拟合一次分布而已。想更快就用groupby.apply效果一样。4.2 不同尺度怎么选SPEI的scale是核心参数不同尺度对应不同类型干旱尺度典型应用生理/水文含义SPEI-1气象干旱、逐月异常当月水分亏缺响应快SPEI-3农业干旱季度水分状况影响作物关键期SPEI-6水资源调度半年累积中小河流流量SPEI-12水文干旱年尺度水分累积影响水库、地下水SPEI-24/48长期干旱演变反映年代际水分背景我常用组合SPEI-3做主图、SPEI-12做趋势背景。两张图放一起能很清晰看出“短期干旱叠加长期持续干旱”的复合风险。4.3 校验与数据源对比如果项目涉及论文或正式报告建议把计算结果和两个数据源对比SPEIbase全球数据集下载对应经纬度格点数据对比同期的SPEI-3曲线国家气候中心/区域气候公报看重大干旱事件年份是否一致。注意数据源不同、PET算法不同、基准期不同都会导致数值有差异。对比时重点是看趋势和干旱等级不是要求完全重合。我遇到过和SPEIbase相关系数只有0.55的站最后查出来是气象站迁移导致气温序列突变这个坑在后面专门讲。5. 可视化把干旱演变画成看懂干旱的图5.1 单站点时间序列图SPEI图我最常用的是“分级着色线图”SPEI值超过0画蓝色低于0画红色低于-1加深低于-1.5再加深。这样一眼就能看出哪几年是重旱。import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(14, 5)) spei df[spei_3].dropna() ax.plot(spei.index, spei.values, lw0.8, colorgray) ax.fill_between(spei.index, 0, spei.values, wherespei 0, color#2a7abf, alpha0.35) ax.fill_between(spei.index, 0, spei.values, where(spei 0) (spei -1), color#d94e28, alpha0.35) ax.fill_between(spei.index, -1, spei.values, where(spei -1) (spei -1.5), color#d94e28, alpha0.7) ax.fill_between(spei.index, -1.5, spei.values, wherespei -1.5, color#a61e17, alpha0.8) ax.axhline(0, colorblack, lw0.8) ax.set_title(SPEI-3 Time Series (Test Station)) plt.tight_layout() plt.show()如果你用Jupyter环境想交互式看具体年月数值可以换成plotlyimport plotly.graph_objects as go fig go.Figure() fig.add_trace(go.Scatter(xspei.index, yspei.values, modelines, nameSPEI-3)) fig.add_hline(y0, line_width1, line_colorblack) fig.add_hrect(y0-2, y1-1.5, fillcolorred, opacity0.2, line_width0) fig.update_layout(titleSPEI-3 Interactive, templateplotly_white) fig.show()5.2 多站点热力图多样品时空演变最直观的是“站点 × 时间”热力图。比如看华北几个站的SPEI-12年际演变# 假设result包含station_id, date, spei字段 pivot result.pivot_table(indexstation_id, columnsresult.index, valuesspei) import matplotlib.pyplot as plt from matplotlib.colors import TwoSlopeNorm norm TwoSlopeNorm(vmin-2, vcenter0, vmax2) fig, ax plt.subplots(figsize(16, 6)) im ax.imshow(pivot.values, aspectauto, cmapRdBu_r, normnorm, extent[pivot.columns.min().year, pivot.columns.max().year, 0, len(pivot.index)]) ax.set_yticks(np.arange(len(pivot.index)) 0.5) ax.set_yticklabels(pivot.index) plt.colorbar(im, axax, labelSPEI) plt.show()这张图特别适合看“区域同步干旱”如果同一列颜色大范围偏红说明这次干旱是全区同时发生的如果只有某一行偏红就是局地干旱。5.3 空间分布简版站点再多一点可以在底图上画气泡图气泡大小表示站点SPEI绝对值颜色红蓝表示干湿。这里不展开复杂插值直接画散点就够了# df_spatial至少包含lon, lat, spei_mean三列 fig, ax plt.subplots(figsize(10, 8)) sc ax.scatter(df_spatial[lon], df_spatial[lat], cdf_spatial[spei_mean], cmapRdBu_r, vmin-2, vmax2, s80, edgecolorsk, linewidths0.5) plt.colorbar(sc, axax, labelSPEI-3 average for 2022 summer) ax.set_xlabel(Longitude) ax.set_ylabel(Latitude) plt.show()如果需要正式地图再叠加cartopy或contextily底图这里给的是最轻量的方案。6. 踩坑记录实际运行中高频出现的问题与解决方案6.1 时间排序与窗口错位这是我自己第一次跑SPEI踩过的坑说出来你可能不信结果图里1985到1995年干旱颜色连成一片最后发现是CSV里年月没有排序rolling窗口把不同季节的数据“混”在一起求和了。rolling是按行顺序滑动的它不管索引的时间语义。所以任何计算之前一定先df df.sort_index()6.2 平年闰年与缺测月pyet.thornthwaite内部处理了闰日但如果你的原始数据是“日值聚合到月值”2月29日会产生0.几毫米的差异。对SPEI这种标准化指数来说这点差异几乎可忽略没必要特殊处理。真正需要注意的是缺测月某月温度缺测PET算出来就是NaND自然也是NaNrolling窗口到那里会少一个数据。我建议缺测比例小于5%的站点用线性插值补齐再算缺测比例超过10%的站点放弃这个站或者只算到某个时期别硬补。6.3 拟合参数异常如果你运行loglogistic_pwm_fit遇到beta 1或者PWM分母接近0先别急着怀疑代码。常见原因有三个序列太短少于30年序列里有大量连0的月份(极端干旱区)分布形态退化数据源有拼接断点某段降水系统性偏差。我遇到过某个站SPEI-12算出来直线下降查到最后是2005年前后站点迁址降水观测值整体变化了。这种数据问题不是统计方法能救的务必回原始资料核实。6.4 可视化中文乱码与刻度过密用matplotlib画中文标题时需要先设置中文字体plt.rcParams[font.sans-serif] [SimHei, Noto Sans CJK SC, Arial Unicode MS] plt.rcParams[axes.unicode_minus] False如果时间跨度长x轴刻度会挤成一团建议用ax.xaxis.set_major_locator(plt.MaxNLocator(10))限制刻度数量或者直接切片只画近20年。最后再分享一个小技巧我每次算完全部站点都会顺手生成一张“SPEI-3标准差检查表”把每个站的SPEI均值和标准差打出来均值偏离0超过0.1、标准差偏离1超过0.2的站一律重点复查。大部分时候是时间索引或者缺测的问题少数时候是站点本身气候态特殊。SPEI计算本身不复杂真正花时间的往往是数据质量和口径验证。你如果刚开始用Python做干旱分析建议先用一个已知站点的数据跑通流程再扩展到批量计算。跑通之后把PET换成Hargreaves或Penman-Monteith也只是替换一行代码的事。数据干净、口径统一后面的分析自然顺手。