如何使用Python计算两个时间序列间随时间变化的时滞

方案选型 #

针对你需要的随时间变化的逐点时滞计算需求，最常用且易实现的是滑动窗口局部互相关法，该方法通过在每个时间点附近取固定窗口，计算窗口内两个序列的互相关峰值对应偏移，得到该点的时滞，完全匹配你给出的时滞定义（相同取值的最小时间差）。

依赖包 #

仅需三个常用Python数据处理包，无需额外安装小众工具：

• pandas：读取处理csv格式数据
• numpy：数组运算
• scipy：互相关计算

完整实现代码 #

1. 导入依赖 #

import pandas as pd
import numpy as np
from scipy.signal import correlate

2. 读取样例数据 #

将你给出的样例数据保存为csv文件后读取：

# 替换为你的数据文件路径
df = pd.read_csv("time_series_data.csv")
time_arr = df["time"].values
series1 = df["series1"].values
series2 = df["series2"].values

3. 时变时滞计算函数 #

def calc_time_varying_lag(time, s_ref, s_align, window_size=21, max_lag_points=10):
    """
    计算两个时间序列的逐点时变时滞
    参数说明：
        time: 时间序列数组，要求采样间隔均匀
        s_ref: 参考序列（以该序列为基准计算时滞）
        s_align: 待对齐序列
        window_size: 滑动窗口大小，建议取奇数，需大于2倍最大时滞点数
        max_lag_points: 允许的最大时滞的点数（单侧），根据数据先验设置
    返回值：
        每个时间点对应的时滞值，单位与time一致：
        正数表示s_align滞后于s_ref，负数表示s_align超前于s_ref
        边界点因窗口不足返回NaN
    """
    n_points = len(time)
    sample_interval = time[1] - time[0]
    lag_result = np.full(n_points, np.nan)
    half_win = window_size // 2

    for center_idx in range(half_win, n_points - half_win):
        # 截取当前窗口的两个序列片段
        win_ref = s_ref[center_idx - half_win : center_idx + half_win + 1]
        win_align = s_align[center_idx - half_win : center_idx + half_win + 1]
        # 标准化片段，消除幅值差异对互相关的影响
        win_ref = (win_ref - win_ref.mean()) / win_ref.std()
        win_align = (win_align - win_align.mean()) / win_align.std()
        # 计算等长互相关
        corr = correlate(win_ref, win_align, mode="same")
        # 仅在允许的时滞范围内找最大相关峰值
        corr_mid = len(corr) // 2
        valid_corr_range = corr[corr_mid - max_lag_points : corr_mid + max_lag_points + 1]
        # 峰值偏移转换为时间单位时滞
        peak_offset = np.argmax(valid_corr_range) - max_lag_points
        lag_result[center_idx] = peak_offset * sample_interval
    return lag_result

4. 调用函数输出结果 #

# 参数可根据你的数据特性调整
lag_arr = calc_time_varying_lag(time_arr, series1, series2, window_size=21, max_lag_points=10)
# 合并时间与时滞结果输出
output_df = pd.DataFrame({
    "时间": time_arr,
    "时滞值": lag_arr
})
# 打印非空结果（去掉边界无法计算的点）
print(output_df.dropna().to_string(index=False))

参数调优说明 #

• 窗口大小：如果你的数据时滞变化很慢，可增大窗口提升结果稳定性；如果时滞变化快，需减小窗口避免平滑掉时变特性
• 最大时滞点数：设置为你预估的最大可能时滞对应的采样点数即可，过大容易引入错误匹配
• 噪声预处理：如果你的序列噪声较大，可先使用Savitzky-Golay滤波或移动平均做平滑处理后再计算时滞，结果会更准确

替代方案 #

如果两个序列存在局部拉伸/压缩的形变，不是单纯的时变平移，可使用滑动窗口动态时间规整（DTW）的方式计算局部对齐偏移，计算复杂度更高但适配性更强。

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