企业官网建设流程全解析

用Python做SPEI（标准化降水蒸散指数）计算，最难的不是公式，而是把数据准备、PET计算、log-logistic拟合和结果可视化串成一条能跑的流水线。尤其当你手头只有一份Excel格式的月降水、月均温表，想直接出一张能写进报告里的SPEI-3曲线图时，网络上的资料往往要么只讲R，要么只贴一段不完整的伪代码。这篇文章针对这个完整链路给出可复现代码和避坑经验，适合气象、水文、农业领域的研究生、工程师和数据分析师。文章里所有代码我都按“能直接复制运行”的标准写，数据格式、参数含义、常见报错也会一并交代清楚。

1. SPEI计算原理与整体设计思路

1.1 SPEI是什么：从“降水单指标”到“气候水平衡”

SPEI全称是Standardized Precipitation Evapotranspiration Index，标准化降水蒸散指数。它由Vicente-Serrano等人在2010年提出，核心思路并不复杂：把每个月的降水减去潜在蒸散（PET），得到一个“水分盈余或亏损”序列，再对这个序列做标准化处理，得到正负值。正值代表湿润，负值代表干旱，数值越小干旱越严重。

相比只考虑降水的SPI，SPEI的最大优势在于把温度变化带了进来。温度升高会直接推高潜在蒸散，导致即使降水不变，水分亏空也可能加剧。这一点在气候变化背景下非常重要。用一个生活化类比来解释：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 matplotlib

pandas负责数据处理和滚动累积，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计算对输入数据有明确要求，最理想的就是一张按月组织的表格，包含三个核心字段：日期、降水、平均气温。

我建议的格式是这样的：

prcp是月降水量，单位必须是毫米；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公式计算潜在蒸散PET

Thornthwaite公式的核心逻辑是：先用多年平均月气温计算一个“热指数I”，再根据当前月气温、月天数和平均日照时数估算该月蒸散能力。公式分三步。

第一步，计算热指数I。假设某站有1991到2020年共30年资料，先求出每个自然月（1月到12月）的多年平均气温 (T_m)，再计算：

[ I = \sum_{m=1}^{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小于等于0，PET直接取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, index=df.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(days=14) 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, index=df.index)

这段代码里有两个要注意的地方。一个是热指数I只用多年平均月气温算一次，不是每个月重算；另一个是N的估算用了天文公式，纬度越高的站点，秋冬月份N越小，PET也会相应降低，更贴近真实蒸散过程。如果不想这么复杂，也可以用平均日照12小时近似，但高纬度地区误差会比较大。

3.2 计算气候水平衡D并完成时间尺度累积

PET拿到以后，气候水平衡就是一行减法：

d = df['prcp'] - pet d.name = 'D'

d的每个值代表该月的“水分盈余”或“水分亏空”。正值说明降水比蒸散多，负值说明入不敷出。

时间尺度累积有两种实现方式，一种是循环累加，一种是pandas的rolling滚动求和。rolling写法更简洁、更不容易错：

scale = 3 spei_running = d.rolling(window=scale).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, dtype=float)) 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，再通过正态分位数函数得到SPEI：

def 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 = f'SPEI-{scale}' return spei

norm.ppf就是把累积概率变成标准正态分布的分位数。比如概率是0.1时，SPEI约等于-1.28，代表发生了累积概率只有10%的干旱事件。

3.4 封装成可直接调用的calculate_spei函数

把前面所有步骤合在一起，封装成一个对使用者友好的函数。以后只要提供降水、气温、纬度和日期索引，就能一次性拿到SPEI序列：

def calculate_spei(prcp, tmean, lat, dates, scale=3): 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( prcp=df['prcp'], tmean=df['tmean'], lat=30.5, dates=df.index, scale=3 ) print(spei_3.tail())

这样整个计算逻辑就闭环了。需要算SPEI-1还是SPEI-12，只需要改scale参数。

4. SPEI结果可视化：时间序列、多尺度对比与正态检验

4.1 SPEI时间序列图与干旱等级阈值

可视化是SPEI分析里最直观的输出，也是写报告时用得最多的部分。最基础的就是画SPEI随时间的折线图，并加上阈值线。

常用干旱等级划分如下：

绘图代码很简单：

fig, ax = plt.subplots(figsize=(12, 5)) spei_3.plot(ax=ax, color='#2c3e50', lw=0.8) ax.axhline(0, color='gray', lw=0.8) ax.axhline(-1.0, color='orange', linestyle='--', lw=0.8) ax.axhline(-1.5, color='red', linestyle='--', lw=0.8) ax.axhline(-2.0, color='darkred', linestyle='--', lw=0.8) ax.fill_between(spei_3.index, -10, 0, color='lightblue', alpha=0.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), sharex=True) for ax, scale in zip(axes, [1, 3, 12]): spei_values = calculate_spei( prcp=df['prcp'], tmean=df['tmean'], lat=30.5, dates=df.index, scale=scale ) spei_values.plot(ax=ax, lw=0.8) ax.axhline(0, color='gray', lw=0.8) ax.axhline(-1, color='orange', linestyle='--', lw=0.6) ax.axhline(-1.5, color='red', linestyle='--', lw=0.6) ax.set_ylabel(f'SPEI-{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, bins=30, density=True, alpha=0.5, label='SPEI-3') ax.plot(x_range, norm.pdf(x_range), color='black', lw=1.2, label='N(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, dist='norm', plot=plt) plt.title('QQ Plot of SPEI-3') plt.show()

QQ图上的点越接近一条直线，说明SPEI序列越接近标准正态分布。一般只要主体部分贴合，两端有些离散是可以接受的。

5. 常见问题与避坑指南

5.1 运行时最常见的四类报错

我在调试SPEI代码过程中，几乎把所有能踩的坑都踩了一遍。下面这四个问题出现的频率最高，而且都比较隐蔽。

其中“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曲线干净漂亮。