SRF腔体LLRF波形数据处理：从三通道校准到动态失谐反演

Posted on 2026-01-05 In note Views: Disqus:

在脉冲超导射频（SRF）腔体的低电平射频（LLRF）系统研究中，如何将采集到的原始波形数据转换为可靠的物理参数是一个核心问题。本文基于一组实际的LLRF测试数据，系统梳理从三通道复数校准、动态失谐提取到品质因数反演的完整流程，并提供可直接运行的Python实现。

问题背景：矢量关系的校准需求

SRF腔体的LLRF系统通常采集三路信号：

前向波（Forward）：forward_phase, forward_amp
腔体场（Cavity）：cavity_phase, cavity_amp
反射波（Reflected）：reflected_phase, reflected_amp

这些测量值与真实物理量之间存在复数增益误差（幅度刻度+相位偏置）。物理上，三路信号必须满足矢量叠加关系：

$$V_c^{\text{true}} = V_f^{\text{true}} + V_r^{\text{true}}$$

我们的目标是通过测量数据标定出各校准系数，使计算值与实测腔体场一致。

三通道全标定方法

数学模型

对每一路测量信号引入独立的复数增益系数：

$$
\begin{cases}
V_f^{\text{true}} = k_f V_f^{\text{raw}}
V_r^{\text{true}} = k_r V_r^{\text{raw}}
V_c^{\text{true}} = k_c V_c^{\text{raw}}
\end{cases}
$$

物理约束条件为：

$$k_c V_c^{\text{raw}} = k_f V_f^{\text{raw}} + k_r V_r^{\text{raw}}$$

为避免尺度不唯一性，工程上采用规范固定：以腔体通道为参考，设k_c = 1。此时问题简化为复数线性最小二乘问题。

标定结果

基于实测数据，复数最小二乘拟合得到：

$$
\begin{aligned}
k_f &= -0.3326 - j0.2327 \quad (\text{Forward通道})
k_r &= 0.1076 - j0.4137 \quad (\text{Reflected通道})
k_c &= 1.0 \quad (\text{Cavity通道})
\end{aligned}
$$

转换为工程常用的幅度×相位形式：

Forward通道系数
$$
|k_f| = 0.405, \quad \angle k_f = -145.0^\circ
$$
校准关系：
$$V_f^{\text{true}} = 0.405 \cdot V_f^{\text{raw}} \cdot e^{-j145.0^\circ}$$

Reflected通道系数
$$
|k_r| = 0.427, \quad \angle k_r = -75.4^\circ
$$

代码实现

# 复数系数构建
deg2rad = np.pi / 180.0
C_f = 0.405 * np.exp(-1j * 145.0 * deg2rad)
C_r = 0.427 * np.exp(-1j * 75.4 * deg2rad)

# 校准应用
Vf_true = C_f * raw_f
Vr_true = C_r * raw_r
Vc_true = raw_c * V_FACTOR  # V_FACTOR为腔体ADC到物理电压转换因子

动态失谐反演

腔体动力学方程

在无束流条件下，腔体包络方程为：

$$\frac{dV_c}{dt} + \left(\omega_{1/2} - j\Delta\omega\right)V_c = 2\omega_{1/2}\frac{\beta}{\beta+1}V_f$$

其中：

$\omega_{1/2} = \omega_0/(2Q_L)$ 为半带宽
$\Delta\omega$ 为动态失谐量
$\beta = Q_0/Q_e$ 为耦合系数

数值求解方法

通过代数变形，可直接计算半带宽和失谐：

def calculate_dynamic_parameters(vc, vf, dt):
    d_vc = np.gradient(vc, dt)
    vc_conj = np.conj(vc)
    vc_abs_sq = np.abs(vc)**2
    
    term_LHS = vc_conj * d_vc
    term_2VfVc = 2 * vf * vc_conj
    
    # 半带宽
    num_w = np.real(term_LHS)
    den_w = np.real(term_2VfVc - vc_abs_sq)
    
    w12 = num_w / den_w
    
    # 失谐量
    dw = (np.imag(term_LHS) - w12 * np.imag(term_2VfVc)) / vc_abs_sq
    
    return w12, dw

关键处理：在场强过低区域（<5%最大值）将计算结果置为NaN，避免数值不稳定。

时间分段标定

在脉冲运行中，LLRF链路的瞬态响应可能导致复数增益随时间变化。推荐采用滑动窗口法：

def segment_calibration(Vf, Vr, Vc, idx):
    A = np.column_stack([Vf[idx], Vr[idx]])
    b = Vc[idx]
    kf, kr = np.linalg.lstsq(A, b, rcond=None)[0]
    return kf, kr

# 滑动窗口计算
window = 100
kf_list, kr_list = [], []
for i in range(len(Vc) - window):
    idx = slice(i, i + window)
    kf, kr = segment_calibration(Vf_raw, Vr_raw, Vc_raw, idx)
    kf_list.append(kf)
    kr_list.append(kr)

通过比较不同时段的系数变化，可判断在线自适应标定的必要性（幅度漂移>2-3%或相位漂移>2-3°时需更新）。

Q0、β反演与能量守恒检查

参数关系

对于单端口SRF腔：

$$
\begin{cases}
\displaystyle \frac{1}{Q_L} = \frac{1}{Q_0} + \frac{1}{Q_e}
\beta = \displaystyle \frac{Q_0}{Q_e}
\end{cases}
$$

可得：
$$Q_L = \frac{Q_0}{1+\beta}, \quad \beta = \frac{Q_0}{Q_L} - 1$$

稳态反射法求β

在平顶区（失谐≈0）使用反射信号：

$$\beta = \frac{1 - |V_r/V_f|}{1 + |V_r/V_f|}$$

# 稳态区取值
mask_steady = (time > 50) & (field > threshold)
Vr_steady = Vr_smooth[mask_steady][-1]
Vf_steady = Vf_smooth[mask_steady][-1]

beta = (1 - np.abs(Vr_steady/Vf_steady)) / (1 + np.abs(Vr_steady/Vf_steady))
Q0 = QL_const * (1 + beta)  # QL_const = 2.78e5

Pickup端口影响评估

当存在弱耦合Pickup端口（$Q_{e,pu}=1.06\times10^{12}$）时：

$$\frac{1}{Q_L} = \frac{1}{Q_0} + \frac{1}{Q_{e,main}} + \frac{1}{Q_{e,pu}}$$

定量判断：
$$\frac{1/Q_{e,pu}}{1/Q_L} \sim \frac{10^{-12}}{10^{-6}} = 10^{-6}$$

Pickup对QL和β的贡献低于$10^{-4}$量级，严格可忽略。这在论文中应明确说明：

The contribution of the pickup probe to the total external loading is negligible due to its extremely weak coupling ($Q_{e,pu}=1.06\times10^{12}$), and was therefore ignored in the dynamic analysis.

能量守恒检查

三通道标定后必须满足：

$$|V_f|^2 = |V_c|^2 \frac{Q_L}{Q_e} + |V_r|^2$$

energy_diff = np.abs(Vf_smooth)**2 - \
              (np.abs(Vc_smooth)**2 * QL_const/Qe_main + np.abs(Vr_smooth)**2)
energy_rel_error = np.abs(energy_diff) / np.abs(Vf_smooth)**2

# 合格标准：最大相对误差 < 1%
print(f"Max relative energy error: {np.nanmax(energy_rel_error):.3e}")

完整处理流程代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter

# ========== 配置参数 ==========
fs = 2e6
dt = 1/fs
f0 = 324e6

# 校准系数（来自三通道标定）
C_f = 0.405 * np.exp(-1j*np.deg2rad(145.0))
C_r = 0.427 * np.exp(-1j*np.deg2rad(75.4))
V_FACTOR = 1038.14

# 腔体参数
QL_const = 2.78e5
Qe_main = 2.78e5
w12_const = 2*np.pi*f0 / (2*QL_const)

# ========== 数据加载与校准 ==========
# raw_f, raw_c, raw_r 为原始复数数据
Vf_true = C_f * raw_f
Vc_true = raw_c * V_FACTOR
Vr_true = C_r * raw_r

# ========== 平滑处理 ==========
window, order = 311, 3
Vc_smooth = savgol_filter(np.real(Vc_true), window, order) + \
            1j * savgol_filter(np.imag(Vc_true), window, order)
Vf_smooth = savgol_filter(np.real(Vf_true), window, order) + \
            1j * savgol_filter(np.imag(Vf_true), window, order)

# ========== 动态失谐计算 ==========
dVc = np.gradient(Vc_smooth, dt)
vc_conj = np.conj(Vc_smooth)
vc_abs_sq = np.abs(Vc_smooth)**2

mask = np.abs(Vc_smooth) > np.max(np.abs(Vc_smooth))*0.05

dw_rad = np.full_like(Vc_smooth, np.nan, dtype=float)
dw_rad[mask] = (np.imag(vc_conj*dVc)[mask] - \
                w12_const * np.imag(2*Vf_smooth*vc_conj)[mask]) / vc_abs_sq[mask]

detuning_Hz = dw_rad / (2*np.pi)

# ========== β与Q0反演 ==========
# 动态计算（需稳态条件）
beta_time = np.full_like(Vc_smooth, np.nan)
beta_time[mask] = (1 - np.abs(Vr_smooth[mask]/Vf_smooth[mask])) / \
                  (1 + np.abs(Vr_smooth[mask]/Vf_smooth[mask]))
Q0_time = QL_const * (1 + beta_time)

# 稳态统计
beta_steady = np.nanmedian(beta_time[mask & (time_axis > 50)])
Q0_steady = np.nanmedian(Q0_time[mask & (time_axis > 50)])

# ========== 能量守恒检查 ==========
energy_rel_error = np.abs(np.abs(Vf_smooth)**2 - 
                          (np.abs(Vc_smooth)**2 * QL_const/Qe_main + 
                           np.abs(Vr_smooth)**2)) / np.abs(Vf_smooth)**2

# ========== 可视化 ==========
fig, axes = plt.subplots(4,1, figsize=(12,14), sharex=True)

axes[0].plot(time_axis, np.abs(Vc_smooth), 'g', label='|Vc|')
axes[0].plot(time_axis, np.abs(Vf_smooth), 'b', label='|Vf|')
axes[0].plot(time_axis, np.abs(Vr_smooth), 'r', label='|Vr|')
axes[0].set_ylabel('Voltage [V]')
axes[0].legend()
axes[0].grid(True)

axes[1].plot(time_axis, detuning_Hz, 'm')
axes[1].set_ylabel('Detuning [Hz]')
axes[1].grid(True)

axes[2].plot(time_axis, beta_time, 'r')
axes[2].set_ylabel(r'$\beta$')
axes[2].grid(True)

axes[3].plot(time_axis, Q0_time, 'b', label='Q0')
axes3_twin = axes[3].twinx()
axes3_twin.plot(time_axis, energy_rel_error, 'k', alpha=0.5, label='Energy error')
axes[3].set_ylabel(r'$Q_0$')
axes3_twin.set_ylabel('Rel Error')
axes[3].set_xlabel('Time [μs]')
axes[3].legend(loc='upper left')
axes3_twin.legend(loc='upper right')
axes[3].grid(True)

plt.tight_layout()
plt.show()

print(f"Steady β: {beta_steady:.3e}, Q0: {Q0_steady:.3e}")
print(f"Max energy error: {np.nanmax(energy_rel_error[mask]):.3e}")

工程实践建议

数据质量优先：平滑参数需要根据信噪比调整，过度平滑会丢失LFD动态信息
阈值选择：5%场强阈值适用于大多数情况，低场区数据应舍弃
稳态区判定：建议自动识别平顶区（如np.diff(|Vc|) < tolerance）
β物理性检查：β应≈1-10（临界耦合到强耦合），若出现β≪1或β≫10需检查标定
能量误差监控：将其作为在线系统的健康检查指标

这套方法已在324MHz DSR腔体上验证，完整覆盖了从原始ADC数据到物理参数的动态反演，可直接用于脉冲SRF系统的实时监测或离线分析。