用Python做SPEI标准化降水蒸散指数计算最难的不是公式而是把数据准备、PET计算、log-logistic拟合和结果可视化串成一条能跑的流水线。尤其当你手头只有一份Excel格式的月降水、月均温表想直接出一张能写进报告里的SPEI-3曲线图时网络上的资料往往要么只讲R要么只贴一段不完整的伪代码。这篇文章针对这个完整链路给出可复现代码和避坑经验适合气象、水文、农业领域的研究生、工程师和数据分析师。文章里所有代码我都按“能直接复制运行”的标准写数据格式、参数含义、常见报错也会一并交代清楚。1. SPEI计算原理与整体设计思路1.1 SPEI是什么从“降水单指标”到“气候水平衡”SPEI全称是Standardized Precipitation Evapotranspiration Index标准化降水蒸散指数。它由Vicente-Serrano等人在2010年提出核心思路并不复杂把每个月的降水减去潜在蒸散PET得到一个“水分盈余或亏损”序列再对这个序列做标准化处理得到正负值。正值代表湿润负值代表干旱数值越小干旱越严重。相比只考虑降水的SPISPEI的最大优势在于把温度变化带了进来。温度升高会直接推高潜在蒸散导致即使降水不变水分亏空也可能加剧。这一点在气候变化背景下非常重要。用一个生活化类比来解释SPI像只盯着你的工资收入SPEI则是同时把每月固定开销也算进去看月末到底能剩下多少钱。降水是收入PET是刚性支出D P - PET就是当月的结余。SPEI的应用场景很广包括气象干旱监测、农业干旱评估、水文水资源规划、森林火险预警等。不同行业关注的时间尺度也不同农业通常关注3个月以内水资源管理则更看重12个月甚至更长的尺度。1.2 标准计算流程拆解SPEI的标准计算流程可以分成四步算PET、算气候水平衡、按时间尺度累积、拟合分布并标准化。第一步是计算每个月的潜在蒸散PET。常用方法有Thornthwaite、Hargreaves、Penman-Monteith等。Thornthwaite方法输入简单只需要月平均气温和纬度适合站点数据Penman-Monteith精度高但需要太阳辐射、风速、湿度等多要素输入数据门槛高。在SPEI最常用的R语言实现中默认就支持Thornthwaite因此这篇教程也以Thornthwaite为主。第二步是计算每个月的水平衡值D公式只有一行[ D_i P_i - PET_i ]这里的P是月降水量PET是同一个月的潜在蒸散量单位都统一为毫米。第三步是时间尺度累积。SPEI可以像SPI一样指定尺度比如SPEI-3就是把连续3个月的D值相加。这一步看似简单但实际操作中经常因为索引错位、缺失值而写错。第四步是对累积后的序列拟合三参数log-logistic分布再通过正态分位数转换把累积概率映射为标准正态值。这一步是SPEI区别于一些简化干旱指数的核心也是最容易出问题的地方。1.3 时间尺度选择为什么重要时间尺度本质上决定了SPEI反映的是“短期土壤墒情”还是“长期水文状态”。SPEI-1反映的是非常短期的水分异常对单个月份的降水、气温波动都很敏感适合监测突发性干旱或农业应急决策但噪声也比较大。SPEI-3是农业干旱最常用的尺度对应一个季度的水分累积能过滤掉单月波动同时又能及时捕捉季节尺度上的干旱信号。SPEI-6和SPEI-12则更多用于水文干旱、地下水补给、水库调度等场景因为这些系统对水分亏缺的响应往往有几个月到一年的滞后。选择尺度时还要考虑数据长度。时间尺度越大有效样本数越少。比如你有30年月值数据共360个点计算SPEI-12时有效累积序列只有349个点虽然足够拟合参数但如果你只有5年数据SPEI-12就只剩49个点拟合结果会非常不稳定。实际项目中我建议至少要有20年以上月值数据少于这个量级时SPEI的统计意义会打折扣。2. Python环境搭建与SPEI输入数据准备2.1 一行命令装好依赖库本次教程的依赖非常轻量核心只有四个库pip install pandas numpy scipy matplotlibpandas负责数据处理和滚动累积numpy做数组计算scipy提供正态分布分位数函数和伽马函数matplotlib负责可视化。Python版本建议3.8以上低于3.8大概率不会报错但pandas新版本对旧版本支持已经越来越弱能升就升。如果你之前没装过这些库安装完以后可以在命令行或Notebook里跑一下import pandas as pd import numpy as np from scipy.stats import norm from scipy.special import gamma as gamma_func import matplotlib.pyplot as plt print(pd.__version__) print(np.__version__)能正常输出版本号说明环境已经通了。2.2 气象数据长什么样才算“合格”SPEI计算对输入数据有明确要求最理想的就是一张按月组织的表格包含三个核心字段日期、降水、平均气温。我建议的格式是这样的dateprcptmean1991-01-0132.5-2.11991-02-0128.30.41991-03-0145.66.8prcp是月降水量单位必须是毫米tmean是月平均气温单位必须是摄氏度。这两个单位不能错错了后面计算结果会完全跑偏。读取数据时有个关键点把date列解析成pandas的DatetimeIndex而且频率要明确是月份的开始。推荐以下写法df pd.read_csv(station_data.csv, parse_dates[date]) df df.set_index(date) df.index pd.DatetimeIndex(df.index) df df.sort_index() # 检查日期是否连续是否存在缺月 print(df.index.min(), df.index.max()) print(df.index[df.index.to_series().diff().dt.days ! 1])最后一行打印出来的是“非连续月份”的位置。如果数据中间缺了某个月后面做滚动求和时会直接出现一大段NaN你不一定第一时间发现。这一步检查值得养成习惯。数据来源方面站点观测可以到国家气象信息中心或国家气候中心申请全球格点再分析数据可以用ERA5月值尺度的变量直接下载即可。无论从哪拿数据都要先确认时间范围和缺失情况。2.3 缺失值、异常值与单位坑数据预处理是整个SPEI计算里最容易被低估的环节。我在实际项目里踩过的坑主要体现在三个方面。第一个坑是缺失值处理。月降水如果只有一两个月缺失我建议用前后月份的线性插值补齐如果连续缺失超过三个月直接补出来的可信度很低宁可舍弃这段序列也不要让插值污染后续的拟合。但要注意如果数据本身只有几年一段三个月的缺失就足以让某个时间尺度上的有效样本显著减少建议重新考虑数据源。第二个坑是异常值。气象观测数据里偶尔会出现极端记录比如某个月降水突然变成平时10倍或者气温出现明显不符合季节变化的数值。SPEI对单个极端值并非完全免疫因为log-logistic拟合会被极端值拉偏。处理异常值时不要简单用“超过多少倍标准差就删除”这种一刀切规则最好结合站点气候背景判断。比如某站夏季月降水极少突然出现一个200mm的记录这未必是错误可能是极端暴雨事件但如果是降水序列里出现负值则几乎可以肯定是数据错误需要处理。第三个坑是单位。ERA5下载的降水一般已经是米需要乘以1000换算成毫米气温默认是开尔文需要减去273.15。不少初学者在数据准备阶段没做单位换算结果PET算出来全是天文数字后面再找问题就要花很长时间。建议在读取数据后立刻统一单位并用describe()检查范围print(df.describe())正常月降水不会为负月均温不会出现上下百度的离谱数值。如果看到明显不符合常识的极值先回看源头数据。3. SPEI核心算法实现PET、水平衡与log-logistic拟合3.1 用Thornthwaite公式计算潜在蒸散PETThornthwaite公式的核心逻辑是先用多年平均月气温计算一个“热指数I”再根据当前月气温、月天数和平均日照时数估算该月蒸散能力。公式分三步。第一步计算热指数I。假设某站有1991到2020年共30年资料先求出每个自然月1月到12月的多年平均气温 (T_m)再计算[ I \sum_{m1}^{12} \left(\frac{\max(T_m, 0)}{5}\right)^{1.514} ]注意气温小于等于0的月份按0处理不参与指数计算。第二步根据热指数I计算指数a[ a 6.75 \times 10^{-7} I^3 - 7.71 \times 10^{-6} I^2 1.792 \times 10^{-2} I 0.49239 ]第三步对每一个具体月份如果当月平均气温T小于等于0PET直接取0如果T大于0则[ PET 16 \times \left(\frac{10 T}{I}\right)^a \times \frac{N}{12} \times \frac{day}{30} ]这里的N是当月平均日照时数单位小时day是当月天数。N可以通过纬度和月中日序用天文学公式估算。完整的Python实现如下def calc_thornthwaite_pet(tmean, lat, dates): # tmean: 月平均气温Series单位°C # lat: 站点纬度单位度 # dates: 与tmean对应的DatetimeIndex df pd.DataFrame({tmean: tmean}) df.index pd.DatetimeIndex(dates) # 1. 多年平均月气温 clim df.groupby(df.index.month)[tmean].mean() clim clim.reindex(range(1, 13)) I_month np.maximum(clim / 5.0, 0) ** 1.514 I I_month.sum() if I 0: return pd.Series(0.0, indexdf.index) # 2. 指数a a (6.75e-7 * I**3 - 7.71e-6 * I**2 1.792e-2 * I 0.49239) # 3. 月中日序与日照时数 mid_month df.index pd.Timedelta(days14) doy mid_month.dayofyear lat_rad np.deg2rad(lat) decl 0.4093 * np.sin(2 * np.pi * (doy - 80) / 365.0) cos_omega -np.tan(lat_rad) * np.tan(decl) cos_omega np.clip(cos_omega, -1, 1) N 24.0 / np.pi * np.arccos(cos_omega) # 月平均日照小时 days_in_month df.index.days_in_month T df[tmean].values pet np.where( T 0, 16.0 * (10.0 * T / I) ** a * (N / 12.0) * (days_in_month / 30.0), 0.0 ) return pd.Series(pet, indexdf.index)这段代码里有两个要注意的地方。一个是热指数I只用多年平均月气温算一次不是每个月重算另一个是N的估算用了天文公式纬度越高的站点秋冬月份N越小PET也会相应降低更贴近真实蒸散过程。如果不想这么复杂也可以用平均日照12小时近似但高纬度地区误差会比较大。3.2 计算气候水平衡D并完成时间尺度累积PET拿到以后气候水平衡就是一行减法d df[prcp] - pet d.name Dd的每个值代表该月的“水分盈余”或“水分亏空”。正值说明降水比蒸散多负值说明入不敷出。时间尺度累积有两种实现方式一种是循环累加一种是pandas的rolling滚动求和。rolling写法更简洁、更不容易错scale 3 spei_running d.rolling(windowscale).sum()这句话的含义是把连续3个月的D相加得到一个累积D序列。前面scale-1个位置是NaN这是正常现象因为样本不够。不同时间尺度用的是同一个d序列只是window参数不同。注意这里不应使用expanding()或cumsum()。cumsum是从序列起点一直累加到最后得不到固定3个月窗口的滑动累积。3.3 三参数log-logistic分布的拟合与标准化这是SPEI算法的核心步骤原理可以这样理解累积后的D序列并不服从正态分布往往带有偏态因此需要先找一个能刻画这种偏态的理论分布去拟合它再把这个分布下的累积概率映射到标准正态分布上得到最终的SPEI值。SPEI原始论文使用三参数log-logistic分布。拟合方法通常有最大似然和L矩两种。L矩方法对干湿序列更稳健也是R语言SPEI包的默认选择所以这里我用L矩实现。L矩方法的核心是先求序列的概率加权矩b0、b1、b2再转成L矩L1、L2、L3最后由L矩解出分布参数。推导过程不展开直接给出封装好的代码from scipy.special import gamma as gamma_func def loglogistic_lmom_fit(data): x np.sort(np.asarray(data, dtypefloat)) n len(x) if n 10: raise ValueError(样本数量太少至少需要10个有效累积值) j np.arange(1, n 1) b0 np.mean(x) b1 np.mean((j - 1) / (n - 1) * x) b2 np.mean((j - 1) * (j - 2) / ((n - 1) * (n - 2)) * x) l1 b0 l2 2 * b1 - b0 l3 6 * b2 - 6 * b1 b0 if l2 0 or l3 0: raise ValueError(L矩无效可能是序列过于均匀或出现极端值) beta l2 / l3 a 1.0 / beta gamma_ratio gamma_func(1 a) * gamma_func(1 - a) alpha l2 / (a * gamma_ratio) gamma_loc l1 - beta * l2 return alpha, beta, gamma_loc代码返回三个参数alpha是尺度参数beta是形状参数gamma是位置参数分别对应log-logistic分布的三个参量。得到参数后对每个累积D值计算CDF再通过正态分位数函数得到SPEIdef loglogistic_cdf(x, alpha, beta, gamma_loc): if x gamma_loc: return 0.0 y (x - gamma_loc) / alpha return 1.0 / (1.0 y**(-beta)) def spei_from_series(series, scale): roll series.rolling(scale).sum().dropna() if len(roll) 10: return roll alpha, beta, gamma_loc loglogistic_lmom_fit(roll) prob roll.apply( lambda x: loglogistic_cdf(x, alpha, beta, gamma_loc) ) prob prob.clip(1e-10, 1 - 1e-10) spei norm.ppf(prob) spei.name fSPEI-{scale} return speinorm.ppf就是把累积概率变成标准正态分布的分位数。比如概率是0.1时SPEI约等于-1.28代表发生了累积概率只有10%的干旱事件。3.4 封装成可直接调用的calculate_spei函数把前面所有步骤合在一起封装成一个对使用者友好的函数。以后只要提供降水、气温、纬度和日期索引就能一次性拿到SPEI序列def calculate_spei(prcp, tmean, lat, dates, scale3): df pd.DataFrame({prcp: prcp, tmean: tmean}) df[date] pd.to_datetime(dates) df df.set_index(date).sort_index() pet calc_thornthwaite_pet(df[tmean], lat, df.index) d df[prcp] - pet spei spei_from_series(d, scale) return spei调用示例spei_3 calculate_spei( prcpdf[prcp], tmeandf[tmean], lat30.5, datesdf.index, scale3 ) print(spei_3.tail())这样整个计算逻辑就闭环了。需要算SPEI-1还是SPEI-12只需要改scale参数。4. SPEI结果可视化时间序列、多尺度对比与正态检验4.1 SPEI时间序列图与干旱等级阈值可视化是SPEI分析里最直观的输出也是写报告时用得最多的部分。最基础的就是画SPEI随时间的折线图并加上阈值线。常用干旱等级划分如下SPEI值区间等级 -2.0极端干旱-2.0 ~ -1.5严重干旱-1.5 ~ -1.0中等干旱-1.0 ~ 1.0接近正常1.0 ~ 1.5中等湿润1.5 ~ 2.0严重湿润 2.0极端湿润绘图代码很简单fig, ax plt.subplots(figsize(12, 5)) spei_3.plot(axax, color#2c3e50, lw0.8) ax.axhline(0, colorgray, lw0.8) ax.axhline(-1.0, colororange, linestyle--, lw0.8) ax.axhline(-1.5, colorred, linestyle--, lw0.8) ax.axhline(-2.0, colordarkred, linestyle--, lw0.8) ax.fill_between(spei_3.index, -10, 0, colorlightblue, alpha0.15) ax.set_title(SPEI-3 Time Series) ax.set_xlabel(Date) ax.set_ylabel(SPEI) ax.legend([SPEI-3]) plt.tight_layout() plt.show()fill_between画的浅蓝色区域能直观区分干湿期。红色阈值线代表干旱等级如果曲线频繁跌破-1.5说明这个区域近年来干旱事件发生频率明显偏高。4.2 多时间尺度对比看短期干旱与长期缺水只画一条SPEI-3曲线可能会遗漏长期累积缺水的信号。更好的做法是把SPEI-1、SPEI-3、SPEI-12放在同一张图里对比让短期波动和长期趋势同时呈现。fig, axes plt.subplots(3, 1, figsize(12, 9), sharexTrue) for ax, scale in zip(axes, [1, 3, 12]): spei_values calculate_spei( prcpdf[prcp], tmeandf[tmean], lat30.5, datesdf.index, scalescale ) spei_values.plot(axax, lw0.8) ax.axhline(0, colorgray, lw0.8) ax.axhline(-1, colororange, linestyle--, lw0.6) ax.axhline(-1.5, colorred, linestyle--, lw0.6) ax.set_ylabel(fSPEI-{scale}) axes[0].set_title(SPEI at Different Time Scales) plt.tight_layout() plt.show()实际看图时可以重点观察SPEI-1往往频繁穿越0轴噪声明显SPEI-12则是一个相对平滑的长周期曲线能够反映出几年的持续湿润或持续干旱。如果短尺度上频繁出现负值事件而长尺度没有明显变化说明干旱是季节性的而非多年级别的缺水。4.3 检验SPEI是否真的接近标准正态分布SPEI在数学上被设计成服从标准正态分布但现实数据拟合后未必完全标准。做一个简单的直方图叠加正态密度可以快速判断结果是否合理。spei_valid spei_3.dropna() x_range np.linspace(-3.5, 3.5, 200) fig, ax plt.subplots(figsize(8, 5)) ax.hist(spei_valid, bins30, densityTrue, alpha0.5, labelSPEI-3) ax.plot(x_range, norm.pdf(x_range), colorblack, lw1.2, labelN(0,1)) ax.set_title(Distribution of SPEI-3) ax.set_xlabel(SPEI) ax.set_ylabel(Density) ax.legend() plt.show()如果直方图明显偏离标准正态曲线比如某个方向拖出很长的尾巴说明数据可能存在极端值或者样本量不足拟合参数可能不够稳健。这时可以再做一次QQ图进一步确认from scipy import stats stats.probplot(spei_valid, distnorm, plotplt) plt.title(QQ Plot of SPEI-3) plt.show()QQ图上的点越接近一条直线说明SPEI序列越接近标准正态分布。一般只要主体部分贴合两端有些离散是可以接受的。5. 常见问题与避坑指南5.1 运行时最常见的四类报错我在调试SPEI代码过程中几乎把所有能踩的坑都踩了一遍。下面这四个问题出现的频率最高而且都比较隐蔽。问题现象常见原因解决办法rolling求和结果全是NaN日期索引不是连续月频或者前面缺月用pd.to_datetime统一索引检查diff()是否每步都是约30天拟合函数报“L矩无效”D序列在某个时间尺度上过于均匀l2或l3非正检查数据是否长时间为常量比如站点长期无降水记录SPEI出现inf或极大负值CDF接近0或1norm.ppf对边界值敏感用clip(1e-10, 1-1e-10)限制累积概率范围PET全部为0热指数I算出来是0或气温全部小于等于0检查气温单位确认不是开尔文确认多年平均气温为正的月份存在其中“L矩无效”最容易让新手懵住。如果某个站点连续12个月降水都为0气温又不高D序列就会非常接近常数。此时log-logistic拟合的L矩会退化SPEI没有统计意义。这种情况下不要强行计算而是应该单独说明该时段数据不适用。5.2 数据细节导致“假干旱”的三个典型场景有时候代码没报错图也画出来了但结果和实际灾情对不上问题往往出在数据细节上。第一个场景是降水量单位看错。ERA5的月降水累计量通常以米为单位如果当成毫米直接用D序列会小到几乎全是负值SPEI会一路跌到-99。这个错误很难靠肉眼发现因为曲线形状看起来还挺正常。解决办法是检查原始变量单位换算后对比一下站点多年平均年降水是否在合理范围。第二个场景是气温数据用了旬值或日平均值而没有聚合成月平均。如果直接把日温当成月温塞进Thornthwaite公式PET会产生巨大波动D序列自然面目全非。正确的做法是先按月聚合df_monthly df_daily.resample(M).agg({ prcp: sum, tmean: mean })第三个场景是站点纬度填错。Thornthwaite公式里的日长修正强烈依赖纬度如果把北纬30度的站写成南纬40度日照时数序列会完全颠倒PET在夏季被严重低估、冬季被严重高估。纬度这个参数写错的时候SPEI基本等于报废。5.3 和R语言SPEI包结果对不上怎么办很多人在算完Python版SPEI后会拿R语言SPEI包的结果做交叉验证发现两者有细微差别就开始怀疑自己写错了。说实话完全一致才奇怪。差异主要来自三个层面。第一是PET计算细节Thornthwaite公式里日照时长估算、月中日序取法、月份天数处理都可能有细微不同第二是log-logistic参数估计方式L矩计算中的权重公式如果实现略有差异参数会有小幅变化第三是标准化时对边界概率的处理方式不同。我的建议是如果两条SPEI曲线的趋势高度吻合干旱事件的开始、结束和强度等级基本一致就可以认为Python实现是可靠的。如果偏差明显优先检查PET计算结果因为PET偏差会直接传导到最终SPEI上。最稳妥的做法是先用一份公开的站点数据同时跑R版和Python版把两张图叠在一起看确认一致后再投入业务使用。这也是我在正式项目里必做的验证步骤。最后再分享一个小技巧计算SPEI之前先把D序列的滚动累积值存一份CSV出来。这个中间文件的价值在于万一后面SPEI结果不对你可以快速定位是PET算错了、累积写错了还是拟合出了问题不用从头再跑一遍。祝各位算出来的SPEI曲线干净漂亮。