3. 基于弱相干脉冲的诱骗态 BB84 参数计算

弱相干脉冲诱骗态 BB84 的增益与误码率模型

弱相干光源（WCS）的光子数服从泊松分布，脉冲平均光子数为 μ 时，含 n 个光子的概率为 P(n)=e^−μμⁿ/n!。多光子脉冲会带来光子数分离（PNS）攻击风险，诱骗态方法通过发送不同强度（信号态 μ、诱骗态 ν、真空态 ω≈0）并分别统计增益与误码率，来估计单光子成分的贡献。本节以 GYS 三强度诱骗态 BB84 为对象，完成增益、误码率与安全密钥率的数值计算。

核心公式与参数约定

设 Alice 发送强度 k∈{μ, ν, ω}，Bob 端探测效率 η_B，暗计数率 p_d（每脉冲每探测器），系统总效率 η=η_B·10^−αL/10，其中 α 为光纤损耗系数（dB/km），L 为传输距离（km）。

• 增益（Yield）：Y_n = 1 − (1−η)ⁿ + p_d（近似，忽略暗计数与信号重合的高阶项）。
• 总增益：Q_k = Σ_n P_k(n)·Y_n = 1 − e^−kη + p_d。
• 单光子误码率：e₁ = (e₀·Y₀ + e_d·η) / Y₁，其中 e₀=0.5 为暗计数误码，e_d 为光学对准误码。
• 总误码率：E_kQ_k = e₀·Y₀ + e_d·(1−e^−kη)。
• 安全密钥率（GLLP 下界）：R ≥ q·[−Q_μ·f(E_μ)·H₂(E_μ) + Q₁^L·(1−H₂(e₁^U))]，q=1/2 为基矢筛选因子，f 为纠错效率，H₂ 为二元熵。

数值计算示例

取 α=0.2 dB/km、η_B=0.1、p_d=1×10⁻⁶、e_d=0.01、μ=0.5、ν=0.1、ω=0，纠错效率 f=1.16。以 L=50 km 为例，η=0.1×10⁻¹=0.01。

• Y₀=p_d=1×10⁻⁶；Y₁=1−(1−0.01)+1×10⁻⁶≈0.010001。
• Q_μ=1−e^−0.5×0.01+1×10⁻⁶≈1−0.995012+0.000001≈0.004989。
• Q_ν=1−e^−0.1×0.01+1×10⁻⁶≈0.0010005。
• E_μQ_μ=0.5×1×10⁻⁶+0.01×(1−e^−0.005)≈5×10⁻⁷+0.01×0.0049875≈5.04×10⁻⁵，故 E_μ≈0.0101。
• e₁=(0.5×1×10⁻⁶+0.01×0.01)/0.010001≈(5×10⁻⁷+1×10⁻⁴)/0.010001≈0.01005。
• Q₁^L 用三强度估计：Q₁^L≈(μ²e^−μ/(μν−ν²))·(Q_νe^ν−Q_ωe^ω−(ν²/μ²)(Q_μe^μ−Q_ωe^ω))，代入得 Q₁^L≈0.0039。
• H₂(0.0101)≈0.0808，H₂(0.01005)≈0.0805。
• R≥0.5×[−0.004989×1.16×0.0808+0.0039×(1−0.0805)]≈0.5×[−0.000468+0.003586]≈0.00156，即约 1.56×10⁻³ 每脉冲。

随 L 增大，η 指数下降，Q_μ 与 Q₁^L 同步衰减，而暗计数项相对占比上升，e₁^U 增大，最终 R 降至 0。上图中曲线与 R=0 的交点即该参数组下的最大安全距离（示意约 130 km，实际值需按上述公式逐点计算）。

可执行计算脚本（Python）

import numpy as np

def h2(x):
    x = np.clip(x, 1e-12, 1-1e-12)
    return -x*np.log2(x) - (1-x)*np.log2(1-x)

def key_rate(L, mu=0.5, nu=0.1, omega=0.0,
             eta_B=0.1, pd=1e-6, ed=0.01, alpha=0.2, f=1.16):
    eta = eta_B * 10**(-alpha*L/10)
    Y0 = pd
    Y1 = 1 - (1-eta) + pd
    Q_mu = 1 - np.exp(-mu*eta) + pd
    Q_nu = 1 - np.exp(-nu*eta) + pd
    Q_om = 2*pd
    E_mu = (0.5*Y0 + ed*(1-np.exp(-mu*eta))) / Q_mu
    # 三强度估计 Q1 下界
    num = Q_nu*np.exp(nu) - Q_om - (nu**2/mu**2)*(Q_mu*np.exp(mu) - Q_om)
    Q1L = mu**2*np.exp(-mu)/(mu*nu - nu**2) * num
    Q1L = max(Q1L, 0.0)
    e1U = (E_mu*Q_mu - 0.5*Y0) / max(Q1L, 1e-12)  # 简化上界
    e1U = min(max(e1U, 0.0), 0.5)
    R = 0.5 * (-Q_mu*f*h2(E_mu) + Q1L*(1 - h2(e1U)))
    return R, Q_mu, E_mu, Q1L, e1U

for L in [0, 25, 50, 75, 100, 125, 150]:
    R, Qm, Em, Q1, e1 = key_rate(L)
    print(f"L={L:3d} km  R={R:.3e}  Q_mu={Qm:.3e}  E_mu={Em:.4f}  Q1L={Q1:.3e}  e1U={e1:.4f}")

运行后观察 R 随 L 单调下降并在某距离处穿越 0，该距离即当前参数下的最大安全传输距离。改变 μ、ν 或 η_B 可观察曲线平移，为后续强度与概率优化提供基线。