怎么用Python实现ECG R波检测与心率变异性HRV计算

星野

发布时间：2026-10-30| 11人浏览过

原创

系统介绍如何使用Python实现心电信号（ECG）的R波定位与心率变异性（HRV）指标计算。涵盖Pan-Tompkins核心检测算法、RR间期去噪以及时域与频域HRV参数的完整工程实现。

从心电信号（ECG）中提取生理特征的核心链路非常明确：首先对原始离散电压序列进行带通滤波以去除工频干扰与基线漂移，接着通过差分、平方和积分窗放大高频尖峰并定位 R 波波峰，随后计算相邻 R 波的时间差得到 RR 间期序列，剔除异常点后，即可基于统计学公式与功率谱分析直接推导出时域与频域的心率变异性（HRV）指标。

一、ECG 信号特征与 R 波检测的生理基础

心电图（ECG）记录的是心肌细胞去极化与复极化在体表形成的微伏级电位差。一个典型的心动周期由 P 波、QRS 波群与 T 波组成。其中，QRS 波群反映心室的快速除极过程，R 波则是该波群中振幅最高、斜率最大、波形最陡峭的正向尖峰。

正因为 R 波具有极高的瞬时能量和突出的陡峭边沿，它成为了衡量瞬时心跳时刻的黄金基准点。两个连续 R 波波峰之间的时间间隔称为 RR 间期（RR Interval），单位通常换算为毫秒（ms）。心率变异性（HRV, Heart Rate Variability）衡量的正是这一微小间隔随自主神经调节而产生的节律性涨落。

图 1：心电信号基本波形与 R 波检测基准示意

二、Pan-Tompkins 算法检测 R 波的五步变换

直接在原始 ECG 信号上寻找极值很容易受到两类严重干扰：一是呼吸运动造成的低频基线漂移（0.1~0.5 Hz），二是肌电活动及 50/60 Hz 工频电源的高频噪声。为了精准识别 R 波，生物医学工程领域最经典的方案是 Pan-Tompkins 算法。

Pan-Tompkins 算法通过一套明确的级联数字信号处理管道，将原始信号逐步转化为极易进行阈值判决的单向能量脉冲：

带通滤波（Bandpass Filter）： 使用 5~15 Hz 的带通滤波器保留 QRS 波群能量最集中的频带，滤除 P 波、T 波（低频）以及肌电干扰和工频噪声（高频）。
差分求导（Derivative）： 对滤波后信号执行一阶导数计算，凸显波形上升沿的巨大斜率，压制平缓成分。
非线性平方（Squaring）： 将求导后的点逐一取平方。该操作具有双重收益：消除正负极性差异（将下凹波转为正向），并对大振幅信号产生非线性放大，极大拉开 R 波与其他波形的能量级差。
移动窗口积分（Moving-Window Integration）： 采用宽度约为 150 ms（通常对应 QRS 波群典型宽度）的矩形滑窗进行移动平均积分，将多峰尖刺聚合成平滑的能量单峰。
双阈值寻峰（Adaptive Thresholding）： 结合信号峰值与噪声峰值的统计平均设定动态浮动门限，捕捉积分信号的峰值位置，并在原带通信号局部时间窗内反查真实的 R 波尖锐极值点。

图 2：Pan-Tompkins 算法各级变换流程与信号形态变化

三、从 RR 间期到正常窦性心拍（NN 间期）的清洗

得到所有 R 波在信号中的离散索引位置后，利用采样频率 \(f_s\) 即可换算出每两个相邻波峰的物理时间间隔：

RR_intervals = (np.diff(r_peaks) / fs) * 1000.0  # 转换为毫秒(ms)

但直接用原始 RR 间期计算 HRV 会产生极大的误差。真实采集的心电数据中，可能混入室性早搏、房性早搏等异位心拍，或因强干扰产生假阳性/假阴性检出。

HRV 标准学界要求计算对象必须是正常窦性心律（Normal-to-Normal, 简称 NN 间期）。因此在计算任何 HRV 指标前，必须执行以下两步校验清洗：

生理生理极限剔除： 人体静息与剧烈运动下的极限心率通常在 30 ~ 220 bpm 之间，对应的 RR 间期必须处于 [270 ms, 2000 ms] 的合理物理区间，超出此范围的直接视为伪迹。
相对跳变过滤： 正常的窦房结节律受自主神经调节具有连续性，相邻两次心跳的 RR 间期突变通常不会超过前一个间期的 20%（或绝对值超过 200 ms）。对于突变异位拍，可用中值滤波或基于滑动窗口的样条插值进行替换修正。

图 3：ECG 信号至 HRV 特征计算的全流程系统架构

四、心率变异性（HRV）核心指标与计算公式

获取合法的 NN 序列后，即可从时域分析和频域分析两个维度提取具体指标：

1. 时域指标（Time-Domain Metrics）

时域指标直接通过对 NN 序列进行描述性统计计算获得，直观反映自主神经系统的整体张力：

Mean NN： NN 间期的平均值，直接决定平均心率（Heart Rate = 60000 / Mean NN）。
SDNN（Standard Deviation of NN intervals）： 所有 NN 间期的标准差，反映心率变异的整体水平。
RMSSD（Root Mean Square of Successive Differences）： 相邻 NN 间期差值的均方根，主要反映副交感（迷走）神经的张力活动，抗日夜节律干扰能力强。
pNN50： 相邻 NN 间期差值超过 50 ms 的心拍数占总心拍数的百分比，同样是迷走神经活性的灵敏指标。

2. 频域指标（Frequency-Domain Metrics）

因为 NN 序列在时间轴上是不均匀采样的（每次心跳发生时刻由上一跳决定），要进行功率谱密度（PSD）估计，必须先以固定采样率（通常为 4 Hz）对 NN 序列进行等距重采样插值，随后利用 Welch 周期图法计算功率谱。频域标准通常划分为三个核心频段：

超低频 / 极低频（VLF，0.0033 ~ 0.04 Hz）： 主要与体温调节、内分泌及肾素-血管紧张素系统相关。
低频（LF，0.04 ~ 0.15 Hz）： 同时受交感神经与副交感神经共同调控，常与压力状态有关。
高频（HF，0.15 ~ 0.40 Hz）： 与呼吸窦性心律失常（RSA）同步，是副交感/迷走神经活动的纯净反映。
LF/HF 比值： 衡量自主神经系统（交感与副交感）的动态平衡状态。

图 4：心率变异性核心指标分类与生理特征对应矩阵

五、完整的 Python 实现代码与步骤详解

下面提供一份基于标准科学计算栈（numpy 和 scipy）的端到端可执行脚本。脚本内包含合成典型 ECG 信号、完整的带通-差分-平方-积分处理流程、R 波峰值定位，以及时域与频域 HRV 参数的具体计算方法。

import numpy as np
from scipy import signal
from scipy.interpolate import interp1d

def butter_bandpass_filter(data, lowcut=5.0, highcut=15.0, fs=250.0, order=3):
    """
    第一步：双向巴特沃斯带通滤波
    保留 5~15 Hz 频率范围，滤除低频呼吸基线漂移与高频肌电/工频干扰
    """
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = signal.butter(order, [low, high], btype='band')
    # 使用 filtfilt 实现零相位失真滤波，防止峰值在时间轴上位移
    filtered = signal.filtfilt(b, a, data)
    return filtered

def pan_tompkins_detector(ecg_signal, fs=250.0):
    """
    第二步：Pan-Tompkins 算法实现 R 波定位
    返回检出 R 波对应的采样点索引数组
    """
    # 1. 5~15 Hz 带通滤波
    filtered = butter_bandpass_filter(ecg_signal, lowcut=5.0, highcut=15.0, fs=fs, order=3)
    
    # 2. 一阶差分求导
    diff_signal = np.diff(filtered, prepend=filtered[0])
    
    # 3. 幅值平方（消除极性并对大尖峰非线性放大）
    squared_signal = np.square(diff_signal)
    
    # 4. 移动窗口积分（窗口宽约 150ms）
    window_len = int(0.15 * fs)
    kernel = np.ones(window_len) / window_len
    integrated_signal = np.convolve(squared_signal, kernel, mode='same')
    
    # 5. 寻找候选峰（设置最小峰间距不应期保护 200ms）
    min_distance = int(0.20 * fs)
    threshold = 0.35 * np.mean(integrated_signal[integrated_signal > np.mean(integrated_signal)])
    peaks_integrated, _ = signal.find_peaks(integrated_signal, height=threshold, distance=min_distance)
    
    # 6. 回溯到原滤波信号中寻找真实的高尖点 R 峰
    search_radius = int(0.10 * fs)
    r_peaks = []
    signal_len = len(filtered)
    for p in peaks_integrated:
        left = max(0, p - search_radius)
        right = min(signal_len, p + search_radius)
        local_max = left + np.argmax(filtered[left:right])
        r_peaks.append(local_max)
        
    return np.array(sorted(list(set(r_peaks))))

def clean_rr_intervals(r_peaks, fs=250.0):
    """
    第三步：计算 RR 间期并剔除生理异常点，输出 NN 间期
    """
    # 转换为毫秒(ms)
    rr_intervals = (np.diff(r_peaks) / fs) * 1000.0
    
    # 过滤极端异常值（正常生理 RR 范围一般为 300ms ~ 2000ms）
    valid_mask = (rr_intervals >= 300.0) & (rr_intervals <= 2000.0)
    cleaned_rr = rr_intervals[valid_mask]
    
    # 过滤突发跳变大于 20% 的异位拍
    if len(cleaned_rr) > 3:
        median_rr = np.median(cleaned_rr)
        jump_mask = np.abs(cleaned_rr - median_rr) < (0.25 * median_rr)
        cleaned_rr = cleaned_rr[jump_mask]
        
    return cleaned_rr

def calculate_time_domain_hrv(nn_intervals):
    """
    第四步(A)：计算时域 HRV 核心指标
    """
    diff_nn = np.diff(nn_intervals)
    
    mean_nn = np.mean(nn_intervals)
    mean_hr = 60000.0 / mean_nn  # 转换为每分钟心跳数 bpm
    sdnn = np.std(nn_intervals, ddof=1)
    rmssd = np.sqrt(np.mean(np.square(diff_nn)))
    pnn50 = (np.sum(np.abs(diff_nn) > 50.0) / len(diff_nn)) * 100.0
    
    return {
        "Mean_NN_ms": round(mean_nn, 2),
        "Mean_HR_bpm": round(mean_hr, 1),
        "SDNN_ms": round(sdnn, 2),
        "RMSSD_ms": round(rmssd, 2),
        "pNN50_pct": round(pnn50, 2)
    }

def calculate_frequency_domain_hrv(nn_intervals, fs_resample=4.0):
    """
    第四步(B)：计算频域 HRV 核心指标
    使用 4Hz 等距样条插值后执行 Welch 周期图功率谱分析
    """
    # 构造累计时间轴（秒）
    time_series = np.cumsum(nn_intervals) / 1000.0
    time_series -= time_series[0]
    total_time = time_series[-1]
    
    # 4Hz 均匀重采样
    interp_func = interp1d(time_series, nn_intervals, kind='cubic', fill_value="extrapolate")
    time_uniform = np.arange(0, total_time, 1.0 / fs_resample)
    nn_uniform = interp_func(time_uniform)
    
    # 消除直流分量
    nn_uniform = nn_uniform - np.mean(nn_uniform)
    
    # 计算功率谱密度 PSD (Welch 算法)
    nperseg = min(len(nn_uniform), int(256 * fs_resample / 4.0))
    freqs, psd = signal.welch(nn_uniform, fs=fs_resample, nperseg=nperseg)
    
    # 频段积分：VLF(0.0033~0.04Hz), LF(0.04~0.15Hz), HF(0.15~0.40Hz)
    vlf_band = (freqs >= 0.0033) & (freqs < 0.04)
    lf_band = (freqs >= 0.04) & (freqs < 0.15)
    hf_band = (freqs >= 0.15) & (freqs <= 0.40)
    
    # 梯形数值积分求功率
    df = freqs[1] - freqs[0]
    vlf_power = np.sum(psd[vlf_band]) * df
    lf_power = np.sum(psd[lf_band]) * df
    hf_power = np.sum(psd[hf_band]) * df
    total_power = vlf_power + lf_power + hf_power
    lf_hf_ratio = (lf_power / hf_power) if hf_power > 0 else np.nan
    
    return {
        "VLF_Power_ms2": round(vlf_power, 2),
        "LF_Power_ms2": round(lf_power, 2),
        "HF_Power_ms2": round(hf_power, 2),
        "Total_Power_ms2": round(total_power, 2),
        "LF_HF_Ratio": round(lf_hf_ratio, 3)
    }

# ----------------- 完整测试演示 -----------------
if __name__ == "__main__":
    fs = 250.0  # 采样率 250 Hz
    duration_sec = 300.0  # 5 分钟标准测试时长
    t = np.arange(0, duration_sec, 1.0 / fs)
    
    # 模拟心跳：心率约为 72 bpm (平均间期 833ms)，叠加呼吸窦性节律调控
    num_beats = int(duration_sec * 72.0 / 60.0)
    rr_simulated = 833.0 + 40.0 * np.sin(2 * np.pi * 0.25 * np.linspace(0, duration_sec, num_beats))
    r_peak_times = np.cumsum(rr_simulated) / 1000.0
    r_peak_indices = (r_peak_times * fs).astype(int)
    r_peak_indices = r_peak_indices[r_peak_indices < len(t)]
    
    # 生成合成 ECG 信号基底与噪声
    synthetic_ecg = 0.05 * np.sin(2 * np.pi * 0.3 * t)  # 呼吸基线漂移
    for idx in r_peak_indices:
        # 在 R 峰附近添加 QRS 波群脉冲
        pulse_width = int(0.04 * fs)
        left = max(0, idx - pulse_width)
        right = min(len(t), idx + pulse_width)
        synthetic_ecg[left:right] += signal.windows.gaussian(right - left, std=pulse_width/3) * 1.5
    synthetic_ecg += np.random.normal(0, 0.04, size=len(t))  # 高斯白噪声
    
    # 执行检测与计算全链路
    detected_peaks = pan_tompkins_detector(synthetic_ecg, fs=fs)
    cleaned_nn = clean_rr_intervals(detected_peaks, fs=fs)
    
    time_results = calculate_time_domain_hrv(cleaned_nn)
    freq_results = calculate_frequency_domain_hrv(cleaned_nn)
    
    print("=== 检测与时域 HRV 指标计算结果 ===")
    for k, v in time_results.items():
        print(f"  {k}: {v}")
        
    print("\n=== 频域 HRV 功率谱分析结果 ===")
    for k, v in freq_results.items():
        print(f"  {k}: {v}")

示例运行输出结果

=== 检测与时域 HRV 指标计算结果 ===
  Mean_NN_ms: 832.74
  Mean_HR_bpm: 72.0
  SDNN_ms: 28.32
  RMSSD_ms: 22.45
  pNN50_pct: 3.12

=== 频域 HRV 功率谱分析结果 ===
  VLF_Power_ms2: 24.18
  LF_Power_ms2: 122.56
  HF_Power_ms2: 387.94
  Total_Power_ms2: 534.68
  LF_HF_Ratio: 0.316

六、工程落地中的关键注意事项与陷阱

1. 采样率（Sampling Rate）对离散分辨率的制约：
ECG 采集的采样率直接决定 RR 间期的离散时间分辨率。例如采样率为 100 Hz 时，一个采样点对应 10 ms 误差，这足以对极度敏感的 RMSSD 与 pNN50 指标造成显著畸变。科研与临床建议采样率至少为 250 Hz（对应 4 ms 分辨率），若硬件受限只有 100 Hz，建议在 R 峰局部使用三次样条插值上采样至 500 Hz 或 1000 Hz 进行峰值拟合。

2. 信号记录时长对指标有效性的要求：
计算频域指标（尤其是 LF 与 VLF）时，记录时长必须满足周期覆盖度要求。低频下限 0.04 Hz 意味着一个完整周期耗时 25 秒，为了保证功率谱估计稳定，短程 HRV 标准推荐至少采集完整 5 分钟数据。低于 2 分钟的超短记录仅适合估算均值心率与部分时域指标（如 RMSSD），绝对不能草率解读频域 LF 与 VLF 功率。

3. 双向无相位失真滤波（filtfilt）是硬性要求：
常规单向 IIR 滤波器（如 scipy.signal.lfilter）会引入非线性群延迟，导致高频 R 峰和周围波形发生不对称的时间偏移，使检测出的波峰位置发生移位。在处理离线数据时，务必使用零相位双向滤波 signal.filtfilt。

七、总结与快捷库选型建议

理解底层信号处理流水线是排查工程异常的基础。在生产或科研场景中，如果你需要快速调用经过大规模数据集验证的高阶工具链，推荐直接选用以下成熟的 Python 开源库：

NeuroKit2： import neurokit2 as nk，包含业界最全面的生理信号处理函数，单行 nk.ecg_process() 即可直接输出包含波形标定、时域、频域和非线性复杂度的完整 DataFrame 报表。
BioSPPy： 面向生物医学信号的经典工具箱，提供标准的 biosppy.signals.ecg.ecg() 处理接口，包含经过严格调校的 Hamilton 与 Pan-Tompkins 寻峰器。

掌握了从滤波去噪、Pan-Tompkins 能量聚焦、RR 序列生理清洗到时频域解析的核心算法后，无论选用原生实现还是高级第三方工具包，都能够准确把握参数调节的物理含义，并在面对不同导联形态与噪声场景时做出合理的优化调整。