21-EFT.WP.Metrology.Sync v1.0 | 第8章 偏移/频偏/抖动估计（鲁棒与序贯）

目录 ／ 文档-技术白皮书 ／ 21-EFT.WP.Metrology.Sync v1.0

第8章 偏移/频偏/抖动估计（鲁棒与序贯）

一句话目标：给出 offset/skew/J 的统一估计与发布口径，结合鲁棒统计与序贯滤波（α-β/Kalman/自适应门控），使在线同步在异常与漂移下仍满足 SLO 并可审计。

I. 范围与对象

1. 适用对象
   • PTP/NTP/SyncE/White Rabbit 等同步会话的端到端 offset、频偏 skew（ppm）与抖动 J 估计。
   • 硬件/软件/混合时间戳，单播/组播/边界时钟场景，批处理与流式在线估计。
2. 输入
   • 事件与时间戳：t1,t2,t3,t4，hwts ∈ {0,1}，会话序号，profile_id。
   • 窗口与时基：tau_mono 上的滑动窗口 W_k，步长 Delta_t。
   • 延迟与非对称模型（来自第7章）：d_up, d_down, asym（用于偏差校正）。
   • 质量侧信号：q_len, rho, J_ref。
3. 输出
   • 估计与不确定度：hat_offset, hat_skew_ppm, hat_J，区间或 U = k * u_c。
   • 伺服注入量：offset_corr（含 asym 修正），过程噪声 Q 与观测噪声 R 的在线自整定。
   • 审计与清单：manifest.sync.est.*（模型、参数、阈值、P95/P99、签名）。

II. 名词与变量

1. 随机量与状态
   • y_k：第 k 次观测的偏移观测值（经 asym 校正）。
   • theta_k：真实偏移；s_k：频偏（d theta / dt）。
   • e_k：测量噪声，w_k：过程噪声；J：抖动（短期偏移散布）。
   • 状态向量 x_k = [ theta_k, s_k ]^T。
2. 时序与窗口
   • t_k：tau_mono 上的时间；Delta_t = t_k - t_{k-1}；窗口 W_k = {y_{k-m+1}..y_k}。
   • 滑动统计：MAD, IQR, P95/P99。
3. 量纲与单位
   unit(theta)="s", unit(J)="s", unit(s)="1"（常以 ppm 发布），dim(*)="[T]" 或无量纲。

III. 公设 P608-*

• P608-1（线性状态—观测）：短时内 theta 一阶可微，theta_{k+1} ≈ theta_k + s_k * Delta_t。
• P608-2（时基一致）：估计在 tau_mono 上进行，对外以 ts 发布并附 offset/skew/J。
• P608-3（鲁棒优先）：先用鲁棒统计得到粗估与尺度，再驱动序贯滤波，异常不进入滤波器累计。
• P608-4（自适应噪声）：Q,R 随运行在线整定，受 J 与残差门控影响。
• P608-5（可追溯与量纲）：check_dim(expr) 必过，估计/门限/预算可复现并签名。
• P608-6（两口径一致承诺）：若观测依赖到达时，必须并行记录 T_arr 两口径与 delta_form。

IV. 最小方程 S608-*

• S608-1（对称口径的偏移与往返）
  offset_sym = ( ( t2 - t1 ) - ( t4 - t3 ) ) / 2
  delay_sym = ( ( t2 - t1 ) + ( t4 - t3 ) ) / 2
  offset_corr = offset_sym + asym
• S608-2（观测与状态演化）
  y_k = theta_k + e_k
  theta_{k+1} = theta_k + s_k * Delta_t + w_{theta,k}
  s_{k+1} = s_k + w_{s,k}
• S608-3（鲁棒尺度与门控）
  MAD = median( | y_i - median(y) | ), sigma_hat ≈ 1.4826 * MAD
  gate: accept y_k iff | y_k - median(W_k) | ≤ k_g * sigma_hat
• S608-4（频偏估计：加权/鲁棒斜率）
  WLS 斜率： s_hat = ( ∑ w_i (t_i - t̄)(y_i - ȳ) ) / ( ∑ w_i (t_i - t̄)^2 )
  Theil–Sen： s_hat = median_{i<j} ( (y_j - y_i) / (t_j - t_i) )
  skew_ppm = 1e6 * s_hat
• S608-5（Huber M-估计）
  设残差 r_i = y_i - ( a + s * t_i )，标准化 u_i = r_i / sigma_hat；
  ψ(u) = u (|u| ≤ c)；ψ(u) = c * sign(u) (|u| > c)；
  迭代解 ∑ ψ(u_i) = 0 与 ∑ t_i ψ(u_i) = 0，更新 a,s。
• S608-6（α-β 滤波器）
  预测： theta_{k|k-1} = theta_{k-1} + s_{k-1} * Delta_t
  残差： r_k = y_k - theta_{k|k-1}
  更新： theta_k = theta_{k|k-1} + α * r_k
  s_k = s_{k-1} + ( β / Delta_t ) * r_k
• S608-7（Kalman 二阶模型）
  A = [[1, Delta_t],[0,1]], H = [1,0]
  预测： x_{k|k-1} = A x_{k-1}, P_{k|k-1} = A P_{k-1} A^T + Q
  增益： K_k = P_{k|k-1} H^T ( H P_{k|k-1} H^T + R )^{-1}
  更新： x_k = x_{k|k-1} + K_k ( y_k - H x_{k|k-1} ), P_k = (I - K_k H) P_{k|k-1}
• S608-8（抖动估计与映射）
  J_k^2 = (1 - beta_J) * J_{k-1}^2 + beta_J * r_k^2
  亦可用差分近似： J_k^2 ≈ 0.5 * E[ ( y_k - y_{k-1} )^2 ]
  与 Allan deviation 的近似映射（短 τ）： sigma_y^2(τ) ≈ ( J^2 ) / ( 2 * τ^2 )
• S608-9（序贯变更检测）
  标准化 z_k = r_k / sigma_hat；
  g_k^+ = max( 0, g_{k-1}^+ + z_k - nu ), g_k^- = max( 0, g_{k-1}^- - z_k - nu )
  触发： max( g_k^+, g_k^- ) ≥ h → 认定漂移/跳变，重置或增大 Q。

V. 统计流程 M60-8（鲁棒与序贯估计闭环）

1. 就绪与校核
   • 对齐 tau_mono，生成 W_k；检查 unit/dim 与时间戳完整性，落 TraceID。
   • 应用第7章的 asym 修正得到观测 y_k。
2. 鲁棒粗估
   • 用 median(W_k) 得到 hat_offset_init；以 MAD 求 sigma_hat 与门控阈值 k_g。
   • 估计 hat_skew_ppm_init（Theil–Sen 或 WLS+Huber）。
3. 初始化序贯滤波
   • 用粗估设置 x_0, P_0；以 J_ref 或 sigma_hat 初始化 R，以历史漂移初始化 Q。
   • 选择 α-β 或 Kalman；硬件时间戳优先（较小 R）。
4. 在线更新与自适应
   • 对新观测执行门控，拒绝异常；更新滤波状态与 J_k。
   • 以残差统计自适应整定 Q,R；当负载 rho 上升时放大 R。
   • 输出 hat_offset/hat_skew_ppm/hat_J 与不确定度 U。
5. 变更检测与处置
   • 运行 CUSUM/GLRT；触发后短暂提高 Q 或重置状态，必要时转入 holdover。
   • 将事件写入审计：change_point, Q/R 调整轨迹。
6. 发布与签名
   • 写入 manifest.sync.est.*：模型、窗口、α/β 或 Q/R、P95/P99、失败计数与签名。
   • 面向伺服输出 offset_corr 与过程参数（见第6章对接）。

VI. 契约与断言

• C60-81（量纲）：check_dim(hat_offset) = [T], check_dim(hat_J) = [T]。
• C60-82（SLO）：|hat_offset|_p99 ≤ tol_offset，|hat_skew_ppm|_p99 ≤ tol_skew_ppm，hat_J_p99 ≤ J_budget。
• C60-83（稳定性）：cond(P_k) ≤ kappa_max，发散时进入降权或重置。
• C60-84（异常率）：outlier_rate ≤ tol_outlier；超阈后切换更强鲁棒（Huber→Tukey biweight）。
• C60-85（两口径一致）：涉及到达时的观测，delta_form_p99 ≤ tol_Tarr。
• C60-86（锁定时间）：T_lock ≤ T_lock_slo，解锁后 T_relock ≤ T_relock_slo。

VII. 实现绑定 I60-*（估计算法与接口）

• prep_observation(exchange, asym_model) -> y_k, Delta_t, meta
• estimate_init_robust(Y, T) -> {hat_offset_init, hat_skew_ppm_init, sigma_hat}
• fit_huber_wls(Y, T, c) -> {a, s, sigma_hat}
• alpha_beta_init(alpha, beta, Q, R, x0, P0) -> filter / kalman_init(A,H,Q,R,x0,P0) -> filter
• filter_step(filter, y_k, Delta_t, gate) -> {x_k, P_k, r_k, J_k}
• adaptive_qr_tuner(residuals, rho, policy) -> {Q', R'}
• detect_change(residuals, nu, h) -> flags
• emit_est_manifest(est, policy) -> manifest.sync.est
• 不变量：reproducible(seed)；门控记录可复演；unit/dim 校核通过；delta_form ≤ tol_Tarr（若适用）。

VIII. 交叉引用

• 时间戳与链路口径：见本卷第4章。
• 协议实现与字段：见本卷第5章。
• 伺服耦合与增益安排：见本卷第6章。
• 延迟与非对称建模：见本卷第7章。
• Allan 家族与频稳口径：见《EFT.WP.Metrology.TimeBase v1.0》第7章。
• 异常与漂移治理：见《Methods.Cleaning v1.0》第8章。
• 不确定度传播：见《Methods.CrossStats v1.0》附录E。

IX. 质量度量与风控

1. 指标
   • offset_p50/p95/p99, skew_ppm_p95/p99, J_p95/p99, outlier_rate, gate_ratio，lock_time, relock_time。
   • 滤波健康：trace(P_k), cond(P_k), R/Q 比，change_events。
2. 风控动作
   • offset_p99 越界 → 降权、提高门控、触发路由/介质复核（联动第7章）。
   • J_p99 飙升 → 增大 R、收紧伺服带宽或启用 WR/TSN 调度。
   • 频偏突变 → 触发 CUSUM，短期增大 Q 并限流；必要时进入 holdover。
   • 异常率高 → 切换更强鲁棒器或延长窗口。

小结