3. 载波相位观测值生成与周跳探测

任务目标

生成一段模拟的GNSS载波相位观测值序列，并在其中人为注入周跳，实现基于高次差分的周跳探测算法，输出探测到的周跳位置与大小。

实现原理

载波相位观测值在无周跳时，其时间序列应保持平滑。对相邻历元做差分（一次差分），若仍有周跳则继续差分。周跳会使差分序列在对应历元处产生显著跳变。本案例采用四阶差分进行探测：对无周跳的静态数据，四阶差分值应在噪声水平附近；当某历元发生周跳时，该历元附近的四阶差分值会出现远超阈值的尖峰。

关键完整代码

import numpy as np

def generate_phase_data(n_epochs=300, dt=1.0, freq=1.57542e9, 
                        phase_noise_std=0.01, cycle_slips=None):
    """
    生成模拟载波相位观测值（单位：周）
    cycle_slips: list of (epoch_index, slip_cycles)
    """
    c = 299792458.0
    # 模拟静态接收机，几何距离恒定
    true_range = 20000000.0  # 20,000 km
    wavelength = c / freq
    
    # 基准相位（周）
    base_phase = true_range / wavelength
    
    # 时间序列
    t = np.arange(n_epochs) * dt
    
    # 模拟钟差漂移引起的缓慢相位变化（斜坡+二次项）
    clock_drift = 1e-9 * t + 2e-12 * t**2  # 秒
    phase_drift = clock_drift * freq  # 转换为周
    
    # 随机噪声（周）
    noise = np.random.normal(0, phase_noise_std, n_epochs)
    
    # 初始相位
    phase = np.full(n_epochs, base_phase) + phase_drift + noise
    
    # 注入周跳
    if cycle_slips:
        for epoch_idx, slip in cycle_slips:
            if epoch_idx < n_epochs:
                phase[epoch_idx:] += slip  # 从该历元起所有后续值都偏移
    
    return phase

def detect_cycle_slips(phase, threshold_factor=8.0):
    """
    基于四阶差分的周跳探测
    返回: [(epoch_index, slip_size), ...]
    """
    # 计算四阶差分
    diff4 = np.diff(phase, n=4)
    
    # 四阶差分的标准差（使用中位数偏差估计鲁棒性更好）
    mad = np.median(np.abs(diff4 - np.median(diff4)))
    sigma = 1.4826 * mad
    
    if sigma < 1e-10:
        sigma = np.std(diff4)
    
    threshold = threshold_factor * sigma
    
    # 找到超过阈值的差分位置
    slip_indices = np.where(np.abs(diff4) > threshold)[0]
    
    # 差分后索引偏移：diff4[i] 对应原始历元 i+4 处的变化
    results = []
    for idx in slip_indices:
        epoch = idx + 4  # 周跳实际发生的历元
        
        # 估计周跳大小：四阶差分值 / 组合系数
        # 四阶差分系数: [1, -4, 6, -4, 1]
        # 若在epoch处发生大小为s的阶跃，则diff4[epoch-4] = s
        slip_size = diff4[idx]
        
        # 避免重复检测同一周跳（相邻历元可能同时超限）
        if not results or epoch - results[-1][0] > 2:
            results.append((epoch, slip_size))
    
    return results

# ===== 运行与验证 =====
if __name__ == "__main__":
    np.random.seed(42)
    
    # 生成数据：在历元100处注入5周周跳，在历元200处注入-3周周跳
    true_slips = [(100, 5.0), (200, -3.0)]
    phase = generate_phase_data(n_epochs=300, cycle_slips=true_slips)
    
    # 探测周跳
    detected = detect_cycle_slips(phase, threshold_factor=8.0)
    
    print("真实周跳:", true_slips)
    print("探测结果:", detected)
    
    # 验证
    correct = 0
    for det_epoch, det_size in detected:
        for true_epoch, true_size in true_slips:
            if abs(det_epoch - true_epoch) <= 1 and abs(det_size - true_size) < 1.0:
                correct += 1
                break
    
    print(f"正确探测: {correct}/{len(true_slips)}")
    if correct == len(true_slips):
        print("测试通过：所有周跳均被正确探测")
    else:
        print("测试未完全通过")

运行方式

python cycle_slip_detection.py

验证方法与判定标准

• 数据生成验证：代码使用固定随机种子（seed=42），保证可复现性。生成300个历元的载波相位数据，噪声标准差为0.01周。
• 周跳注入：在第100历元注入+5周跳变，第200历元注入-3周跳变。注入方式为从该历元起所有后续相位值整体偏移，模拟真实周跳的持续性影响。
• 探测算法：四阶差分将周跳转换为单历元尖峰。阈值设为8倍MAD估计的标准差，兼顾灵敏度与抗噪声能力。
• 通过条件：程序输出中“正确探测”数量等于2，且显示“测试通过”。探测到的周跳历元与真实值偏差不超过1个历元，大小偏差小于1周。

关键参数说明

• 四阶差分：对静态或低动态场景，四阶差分可有效消除多项式趋势，使周跳信号凸显。
• 阈值因子：8倍σ在噪声标准差0.01周时，虚警率极低；若噪声增大，可适当降低阈值因子。
• 周跳大小估计：四阶差分值直接等于周跳大小（单位：周），无需额外换算。