如何使用Python pandas库计算时间序列月滞后互相关

核心逻辑：先将日尺度数据按自然月聚合为月尺度序列，再通过序列移位对齐的方式逐滞后步长计算皮尔逊相关系数，全程仅用Pandas内置方法即可实现，无需额外依赖专业信号处理库。

步骤1：日尺度数据转月尺度 #

两个变量物理意义不同，聚合规则不能搞混：

• 降雨量是时段累计通量，按月做求和统计
• 地下水位（GWL）是瞬时观测的状态量，按月做算术平均统计
  直接用Pandas原生resample('M')实现自然月聚合，不要手动按30天滑窗统计，避免跨月对齐误差。

步骤2：序列预处理消除伪相关 #

别跳过这步直接算原始序列的相关：两个序列都存在雨季偏高、旱季偏低的固定季节周期，还有可能存在长期上升/下降趋势，会算出没有实际物理意义的虚高相关值，无法反映真实的入渗滞后规律。
通用预处理流程：

1. 计算月距平：每个月的数值减去该月份的多年平均值，消除季节周期和长期趋势干扰
2. 对距平序列做Z-score标准化，保证最终算出的互相关系数等价于皮尔逊相关系数，取值固定在[-1,1]区间，可直接横向对比大小。

步骤3：逐滞后步长计算互相关 #

明确定义滞后规则：lag=k代表当月GWL与k个月前的降雨量的相关系数，也就是降雨发生k个月后，和GWL的相关程度。
通过Pandas的shift(k)方法将降雨序列向后移动k个位置，和GWL序列对齐后调用内置.corr()方法，即可得到对应滞后步长的相关系数。一般设置最大计算滞后为12个月即可覆盖绝大多数浅层地下水的响应时长，深层地下水可适当调大到24~36个月。

步骤4：识别最优滞后月份 #

所有滞后步长的计算结果存入DataFrame后，取相关系数最大值对应的滞后值，就是GWL对降雨响应的高相关滞后月份。

完整可运行代码 #

import pandas as pd
import numpy as np

# 1. 加载数据，替换为你本地的文件路径
# 要求原始数据包含三列：date(日期格式), rainfall(日降雨量), gwl(日地下水位观测值)
df = pd.read_csv('your_rain_gwl_data.csv', parse_dates=['date'], index_col='date')

# 2. 日尺度转月尺度
monthly_df = pd.DataFrame()
monthly_df['rainfall'] = df['rainfall'].resample('M').sum()  # 月累计降雨
monthly_df['gwl'] = df['gwl'].resample('M').mean()           # 月平均地下水位
monthly_df = monthly_df.dropna()  # 剔除存在缺测的月份

# 3. 序列预处理
# 计算月距平
monthly_df['rain_anom'] = monthly_df.groupby(monthly_df.index.month)['rainfall'].transform(
    lambda x: x - x.mean()
)
monthly_df['gwl_anom'] = monthly_df.groupby(monthly_df.index.month)['gwl'].transform(
    lambda x: x - x.mean()
)
# Z-score标准化
rain_norm = (monthly_df['rain_anom'] - monthly_df['rain_anom'].mean()) / monthly_df['rain_anom'].std()
gwl_norm = (monthly_df['gwl_anom'] - monthly_df['gwl_anom'].mean()) / monthly_df['gwl_anom'].std()

# 4. 逐滞后步长计算互相关
max_lag = 12  # 可根据研究区地下水埋深调整
lag_results = []
for lag in range(0, max_lag + 1):
    # 降雨序列后移lag位，和当前GWL对齐计算相关
    corr_val = rain_norm.shift(lag).corr(gwl_norm)
    lag_results.append({'lag_month': lag, 'corr_coef': corr_val})

corr_table = pd.DataFrame(lag_results)

# 5. 提取最高相关对应的滞后月份
best_match = corr_table.loc[corr_table['corr_coef'].idxmax()]
print(f"最高互相关系数为{best_match['corr_coef']:.3f}，对应GWL对降雨的响应滞后为{int(best_match['lag_month'])}个月")

常见踩坑提示 #

• 移位方向不要写反：如果对GWL序列做shift再和降雨算相关，得到的是GWL超前于降雨的相关，不符合降雨入渗补给的因果逻辑。
• 原始日数据如果存在连续缺测，先做线性插值或者直接剔除缺测时段再聚合，避免移位对齐后出现大量空值干扰计算结果。
• 如果研究区存在强烈的地下水开采、灌溉等人为干扰，可以在预处理阶段增加线性去趋势步骤，进一步剔除人为活动导致的GWL趋势影响。

内容的提问来源于stack exchange，提问作者Irum Khan