V2I 信道模型教学演示 - 3GPP TR 38.901 路径损耗 + Nakagami-m 衰落 + Doppler + 链路持续 (独立脚本，无工程依赖)

v2i-channel-model-demo.py

#!/usr/bin/env python3
"""
V2I 信道模型教学演示脚本
========================

本脚本完整实现了 V2I (Vehicle-to-Infrastructure) 通信系统的信道模型，
适合新人理解和学习 5G V2X 无线通信基础。

信道模型基于:
  - 3GPP TR 38.901 v16.1.0 (路径损耗)
  - Nakagami-m 衰落分布 (小尺度衰落)
  - Doppler 效应 (高速移动惩罚)
  - 链路持续时间建模 (覆盖连续性)

依赖: pip install numpy scipy matplotlib
运行: python v2i-channel-model-demo.py

Author: JLU-MCNS-MEC
"""

import numpy as np
from scipy.stats import nakagami
from scipy.special import gammainc
import matplotlib.pyplot as plt
import matplotlib

matplotlib.rcParams["font.family"] = ["WenQuanYi Micro Hei", "SimHei", "sans-serif"]
matplotlib.rcParams["axes.unicode_minus"] = False


# ============================================================
# 第一部分: 通信系统参数
# ============================================================
# 这些参数定义了 V2I 无线信道的物理特性

# --- 发射端 (车辆) ---
P_tx = 0.3                  # 发射功率 (W)，约 25 dBm
h_veh = 1.5                 # 车辆天线高度 (m)

# --- 接收端 (RSU/MEC) ---
h_rsu = 10.0                # RSU 天线高度 (m)
R_coverage = 150.0          # 覆盖半径 (m)

# --- 无线信道 ---
f_c = 5.9e9                 # 载波频率 (Hz)，5.9 GHz 是 C-V2X 标准频段
c = 3e8                     # 光速 (m/s)
BW = 5e7                    # 总带宽 (Hz)，50 MHz
N_0 = 0.8e-13               # AWGN 噪声功率 (W)
sigma_sf = 6.0              # 阴影衰落标准差 (dB)

# --- Nakagami-m 衰落参数 ---
# m 越大表示衰落越轻微 (m=1 等价于 Rayleigh，m→∞ 等价于无衰落)
m_LoS = 4                   # 视距 (LoS) 信道 shape 参数
m_NLoS = 2                  # 非视距 (NLoS) 信道 shape 参数
nakagami_avg_gain = 1       # 平均功率增益

# --- Doppler ---
doppler_factor = 100.0      # Doppler 惩罚系数

# --- SNR 阈值 ---
snr_threshold_dB = 10.0     # 中断 SNR 阈值 (dB)
snr_threshold = 10 ** (snr_threshold_dB / 10)  # 线性值


# ============================================================
# 第二部分: 信道模型核心函数
# ============================================================

def path_loss_dB(d):
    """
    3GPP TR 38.901 路径损耗模型 (dB)

    分为近场和远场两个区间:
    - 近场 (d ≤ d_bp): 自由空间传播，损耗随 log₁₀(d) 线性增长
    - 远场 (d > d_bp):  多径反射加剧，损耗更快增长

    断点距离 d_bp 取决于天线高度和载波频率。

    参考: 3GPP TR 38.901 Table 7.4.1-1 (UMi - Street Canyon LoS)
    """
    d = np.maximum(d, 0.1)  # 避免 log(0)

    # 断点距离: 一阶菲涅尔区与地面反射交汇点
    d_bp = 4 * h_rsu * h_veh * (f_c / 1e9) / (c / 1e9)

    # 近场: 自由空间损耗为主
    pl_near = 32.4 + 21 * np.log10(d) + 20 * np.log10(f_c / 1e9) + sigma_sf

    # 远场: 多径效应增强，衰减更快
    pl_far = (32.4 + 40 * np.log10(d) + 20 * np.log10(f_c / 1e9)
              - 9.5 * np.log10(d**2 + (h_rsu - h_veh)**2) + sigma_sf)

    return np.where(d <= d_bp, pl_near, pl_far)


def path_loss_nlos_dB(d):
    """
    NLoS 路径损耗 (dB)

    非视距场景下损耗更大，取 LoS 路径损耗和 NLoS 经验公式的较大值。
    """
    pl_los = path_loss_dB(d)
    d = np.maximum(d, 0.1)
    pl_nlos = 35.3 * np.log10(d) + 22.4 + 21.3 * np.log10(f_c / 1e9) - 0.3 * (h_veh - 1.5)
    return np.maximum(pl_los, pl_nlos)


def prob_los(d):
    """
    视距 (LoS) 概率

    距离越远，被建筑物遮挡的可能性越大，LoS 概率越低。
    d ≤ 18m 时假设必然 LoS。
    """
    d = np.asarray(d, dtype=float)
    p = np.where(
        d <= 18,
        1.0,
        18.0 / d + np.exp(-d / 36.0) * (1.0 - 18.0 / d)
    )
    return p


def compute_snr(d):
    """
    信噪比 (SNR) 计算

    综合 LoS 和 NLoS 路径，按各自概率加权:
      channel_gain = P_LoS × G_LoS + (1 - P_LoS) × G_NLoS
      SNR = P_tx × channel_gain / N_0

    其中 G 包含小尺度衰落 (Nakagami) 和大尺度路径损耗。
    """
    d = np.maximum(d, 1e-6)

    # 小尺度衰落: Nakagami-m 分布采样
    small_fading = nakagami(m_LoS, loc=0, scale=np.sqrt(nakagami_avg_gain))
    g_small = np.mean(small_fading.rvs(size=1) ** 2)

    # 大尺度路径损耗 (dB → 线性)
    g_los_large = 10 ** (-path_loss_dB(d) / 10)
    g_nlos_large = 10 ** (-path_loss_nlos_dB(d) / 10)

    # 综合信道增益
    p = prob_los(d)
    g_total = p * (g_small * g_los_large) + (1 - p) * (g_small * g_nlos_large)

    return P_tx * g_total / N_0


def outage_probability(d):
    """
    中断概率

    当 SNR 低于阈值时通信中断。使用 Nakagami-m 分布的 CDF:
      P_out = P_LoS × P(SNR_LoS < γ) + (1-P_LoS) × P(SNR_NLoS < γ)

    其中 P(SNR < γ) = γ_inc(m, m × γ / SNR_avg)，γ_inc 是正则化下不完全 Gamma 函数。
    """
    d = np.maximum(d, 1e-6)

    # LoS / NLoS 的平均 SNR
    g_los = 10 ** (-path_loss_dB(d) / 10)
    g_nlos = 10 ** (-path_loss_nlos_dB(d) / 10)
    snr_los = P_tx * g_los / N_0
    snr_nlos = P_tx * g_nlos / N_0

    # Nakagami-m 中断概率 (正则化下不完全 Gamma 函数)
    x_los = m_LoS * snr_threshold / np.maximum(snr_los, 1e-20)
    p_out_los = gammainc(m_LoS, x_los)

    x_nlos = m_NLoS * snr_threshold / np.maximum(snr_nlos, 1e-20)
    p_out_nlos = gammainc(m_NLoS, x_nlos)

    p = prob_los(d)
    return p * p_out_los + (1 - p) * p_out_nlos


def doppler_penalty(v):
    """
    Doppler 效应惩罚因子

    车辆移动导致载波频率偏移 (多普勒频移):
      f_d = v × f_c / c

    频移越大，信道相干时间越短，有效传输速率下降。
    惩罚模型: α = exp(-factor × 2f_d / BW)
    """
    v = np.asarray(v, dtype=float)
    f_d = v * f_c / c                    # 多普勒频移 (Hz)
    B_d = 2 * f_d                        # 多普勒扩展 (Hz)
    alpha = np.exp(-doppler_factor * B_d / BW)
    return np.clip(alpha, 0.0, 1.0)


def link_duration(d, v):
    """
    链路持续时间

    车辆在 MEC 覆盖区内的停留时间:
      T_link = 2 × √(R² - d²) / v

    假设车辆沿弦方向穿越圆形覆盖区。
    """
    d = np.asarray(d, dtype=float)
    v = np.maximum(np.asarray(v, dtype=float), 0.1)
    inside = d < R_coverage
    remaining = np.where(inside, 2 * np.sqrt(np.maximum(R_coverage**2 - d**2, 0)), 0)
    t = np.clip(remaining / v, 0, 300)
    return t


def link_duration_penalty(d, v):
    """
    链路持续惩罚

    连接时间太短意味着任务可能无法在断链前完成:
      < 5s  → 指数衰减惩罚 (严重)
      5-20s → 线性插值 0.5~1.0
      20-60s → Sigmoid 过渡 0.9~1.0
      ≥ 60s → 无惩罚
    """
    t = link_duration(d, v)
    t = np.asarray(t, dtype=float)

    penalty = np.where(
        t <= 0, 0.0,
        np.where(t < 5, 1.0 - np.exp(-t / 2.0),
        np.where(t < 20, 0.5 + 0.5 * (t - 5) / 15.0,
        np.where(t < 60, 0.9 + 0.1 / (1 + np.exp(-0.1 * (t - 40))),
        1.0))))
    return np.clip(penalty, 0, 1)


def transmission_rate(d, v=0.0, bw_ratio=0.5, n_vehicles=1):
    """
    有效传输速率 (核心公式)

    综合考虑所有信道效应:
      Rate = B_eff × log₂(1 + SNR) × (1 - P_out) × α_doppler × β_link

    其中:
      B_eff = BW × bw_ratio / n_vehicles  (公平带宽分配)
    """
    B_eff = BW * bw_ratio / max(n_vehicles, 1)
    snr = np.maximum(compute_snr(d), 1e-6)
    p_out = outage_probability(d)
    dp = doppler_penalty(v)
    lp = link_duration_penalty(d, v)

    base_rate = B_eff * np.log2(1 + snr)
    return base_rate * (1 - p_out) * dp * lp


def transmission_delay(d, task_bits, v=0.0, bw_ratio=0.5):
    """传输时延 = 任务大小 / 传输速率"""
    rate = transmission_rate(d, v, bw_ratio)
    return np.where(rate > 0, task_bits / rate, 1e10)


# ============================================================
# 第三部分: 可视化演示
# ============================================================

def demo_all():
    """运行全部演示，生成 6 张子图"""
    distances = np.linspace(1, 300, 300)
    speeds = np.linspace(0, 30, 100)

    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    fig.suptitle("V2I 信道模型教学演示 (3GPP TR 38.901)", fontsize=14, fontweight="bold")

    # ---- 图1: 路径损耗 vs 距离 ----
    ax = axes[0, 0]
    ax.plot(distances, path_loss_dB(distances), "b-", label="LoS", linewidth=2)
    ax.plot(distances, path_loss_nlos_dB(distances), "r--", label="NLoS", linewidth=2)
    d_bp = 4 * h_rsu * h_veh * (f_c / 1e9) / (c / 1e9)
    ax.axvline(d_bp, color="gray", linestyle=":", alpha=0.7, label=f"断点 d_bp={d_bp:.0f}m")
    ax.set_xlabel("距离 (m)")
    ax.set_ylabel("路径损耗 (dB)")
    ax.set_title("① 路径损耗 vs 距离")
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

    # ---- 图2: LoS 概率 + 中断概率 ----
    ax = axes[0, 1]
    ax2 = ax.twinx()
    l1, = ax.plot(distances, prob_los(distances), "g-", linewidth=2, label="P(LoS)")
    l2, = ax2.plot(distances, outage_probability(distances), "r-", linewidth=2, label="P(中断)")
    ax.set_xlabel("距离 (m)")
    ax.set_ylabel("LoS 概率", color="g")
    ax2.set_ylabel("中断概率", color="r")
    ax.set_title("② LoS 概率 & 中断概率")
    ax.legend(handles=[l1, l2], fontsize=9)
    ax.grid(True, alpha=0.3)

    # ---- 图3: SNR vs 距离 ----
    ax = axes[0, 2]
    snr_vals = compute_snr(distances)
    snr_dB = 10 * np.log10(np.maximum(snr_vals, 1e-20))
    ax.plot(distances, snr_dB, "b-", linewidth=2)
    ax.axhline(snr_threshold_dB, color="r", linestyle="--", label=f"中断阈值 {snr_threshold_dB} dB")
    ax.set_xlabel("距离 (m)")
    ax.set_ylabel("SNR (dB)")
    ax.set_title("③ 信噪比 vs 距离")
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

    # ---- 图4: Doppler 惩罚 vs 速度 ----
    ax = axes[1, 0]
    ax.plot(speeds, doppler_penalty(speeds), "m-", linewidth=2)
    typical_speeds = [8.3, 13.9, 22.2]  # 30, 50, 80 km/h
    for vs in typical_speeds:
        p = doppler_penalty(vs)
        ax.annotate(f"{vs*3.6:.0f} km/h\n({p:.3f})", xy=(vs, p),
                    fontsize=8, ha="center", va="top",
                    arrowprops=dict(arrowstyle="->", color="gray"))
    ax.set_xlabel("速度 (m/s)")
    ax.set_ylabel("惩罚因子")
    ax.set_title("④ Doppler 惩罚 vs 车速")
    ax.grid(True, alpha=0.3)

    # ---- 图5: 传输速率 (不同速度) ----
    ax = axes[1, 1]
    for v, color, ls in [(0, "b", "-"), (10, "g", "--"), (20, "r", ":")]:
        rates = transmission_rate(distances, v=v, bw_ratio=0.5) / 1e6
        label = f"v={v} m/s ({v*3.6:.0f} km/h)"
        ax.plot(distances, rates, color=color, linestyle=ls, linewidth=2, label=label)
    ax.axvline(R_coverage, color="gray", linestyle=":", alpha=0.5, label=f"覆盖边界 {R_coverage}m")
    ax.set_xlabel("距离 (m)")
    ax.set_ylabel("传输速率 (Mbps)")
    ax.set_title("⑤ 有效传输速率 vs 距离")
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

    # ---- 图6: 链路持续时间 & 惩罚 ----
    ax = axes[1, 2]
    ax2 = ax.twinx()
    for v, color in [(5, "b"), (10, "g"), (20, "r")]:
        t_link = link_duration(distances, v)
        penalty = link_duration_penalty(distances, v)
        ax.plot(distances, t_link, color=color, linestyle="-", linewidth=1.5,
                label=f"持续 v={v}m/s")
        ax2.plot(distances, penalty, color=color, linestyle="--", linewidth=1.5, alpha=0.6)
    ax.set_xlabel("距离 (m)")
    ax.set_ylabel("链路持续 (s)")
    ax2.set_ylabel("惩罚因子 (虚线)")
    ax.set_title("⑥ 链路持续时间 & 惩罚")
    ax.legend(fontsize=8, loc="upper right")
    ax.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig("v2i_channel_model_demo.png", dpi=150, bbox_inches="tight")
    print("图表已保存: v2i_channel_model_demo.png")
    plt.show()


# ============================================================
# 第四部分: 数值验证
# ============================================================

def verify_parameters():
    """打印关键参数并验证典型场景下的通信性能"""
    print("=" * 60)
    print("V2I 信道模型参数总览")
    print("=" * 60)

    print(f"\n【发射端】")
    print(f"  发射功率:     {P_tx} W ({10*np.log10(P_tx*1000):.1f} dBm)")
    print(f"  车辆天线高度: {h_veh} m")

    print(f"\n【接收端】")
    print(f"  RSU 天线高度: {h_rsu} m")
    print(f"  覆盖半径:     {R_coverage} m")

    print(f"\n【信道特性】")
    print(f"  载波频率:     {f_c/1e9} GHz (C-V2X)")
    print(f"  系统带宽:     {BW/1e6} MHz")
    print(f"  噪声功率:     {N_0:.1e} W")
    print(f"  阴影衰落:     {sigma_sf} dB")

    print(f"\n【衰落模型】 Nakagami-m")
    print(f"  LoS  m = {m_LoS}  (衰落轻微，接近 Rice)")
    print(f"  NLoS m = {m_NLoS}  (衰落较强，接近 Rayleigh)")

    d_bp = 4 * h_rsu * h_veh * (f_c / 1e9) / (c / 1e9)
    print(f"\n【路径损耗】 3GPP TR 38.901")
    print(f"  断点距离:     {d_bp:.1f} m")
    print(f"  近场指数:     2.1 (∝ d^2.1)")
    print(f"  远场指数:     4.0 (∝ d^4.0)")

    print(f"\n{'=' * 60}")
    print(f"典型场景验证")
    print(f"{'=' * 60}")

    print(f"\n{'距离':>6s} {'PL(dB)':>8s} {'SNR(dB)':>8s} {'P_out':>8s} "
          f"{'Rate@0':>10s} {'Rate@15':>10s} {'T_link@10':>10s}")
    print("-" * 70)

    for d in [10, 30, 50, 80, 100, 130, 150, 200]:
        pl = path_loss_dB(d)
        snr = 10 * np.log10(max(compute_snr(d), 1e-20))
        p_out = outage_probability(d)
        r0 = transmission_rate(d, v=0, bw_ratio=0.5) / 1e6
        r15 = transmission_rate(d, v=15, bw_ratio=0.5) / 1e6
        t = link_duration(d, 10)
        print(f"{d:>5d}m {pl:>8.1f} {snr:>8.1f} {p_out:>8.4f} "
              f"{r0:>8.2f} Mb {r15:>8.2f} Mb {t:>8.1f}s")

    # 传输时延对比
    task_mb = 0.5
    task_bits = task_mb * 8 * 1e6
    print(f"\n任务大小: {task_mb} MB ({task_bits/1e6:.0f} Mbit)")
    print(f"\n{'距离':>6s} {'时延@v=0':>12s} {'时延@v=10':>12s} {'时延@v=20':>12s}")
    print("-" * 50)
    for d in [20, 50, 100, 140]:
        d0 = transmission_delay(d, task_bits, v=0, bw_ratio=0.5) * 1000
        d10 = transmission_delay(d, task_bits, v=10, bw_ratio=0.5) * 1000
        d20 = transmission_delay(d, task_bits, v=20, bw_ratio=0.5) * 1000
        print(f"{d:>5d}m {d0:>10.2f} ms {d10:>10.2f} ms {d20:>10.2f} ms")

    # 带宽共享影响
    print(f"\n带宽共享影响 (d=80m, v=10m/s, bw_ratio=0.5)")
    print(f"{'车辆数':>6s} {'速率':>12s} {'时延':>12s}")
    print("-" * 35)
    for n in [1, 3, 5, 10]:
        r = transmission_rate(80, v=10, bw_ratio=0.5, n_vehicles=n) / 1e6
        t = task_bits / (r * 1e6) * 1000 if r > 0 else float("inf")
        print(f"{n:>5d}  {r:>10.2f} Mb {t:>10.2f} ms")


# ============================================================
# 主程序
# ============================================================

if __name__ == "__main__":
    verify_parameters()
    print("\n正在生成可视化图表...")
    demo_all()