主题070：城市与室内电波传播

一、引言

电波传播是无线通信系统设计的核心问题之一。随着5G/6G移动通信、物联网(IoT)、智能交通系统(ITS)等技术的快速发展，对城市和室内环境中电波传播特性的精确建模和预测变得愈发重要。城市峡谷中的多径效应、建筑物的穿透损耗、室内复杂的结构布局等因素都会显著影响信号覆盖质量和通信性能。

本教程将系统介绍城市与室内电波传播的基本理论、建模方法和仿真技术，包括射线追踪算法、主导路径模型、统计信道模型等内容，并通过Python代码实现多种传播场景的仿真分析。

二、电波传播基础理论

2.1 自由空间传播模型

自由空间传播是最基本的电波传播模型，假设发射机和接收机之间没有任何障碍物。接收功率可由Friis传输方程描述：

$$P_r=P_t{G}_t{G}_r{\left(\frac{\lambda }{4\pi d}\right)}^2$$

其中：

• $P_t$：发射功率(W)
• $P_r$：接收功率(W)
• $G_t$：发射天线增益
• $G_r$：接收天线增益
• $\lambda$：波长(m)
• $d$：传播距离(m)

自由空间路径损耗(以dB为单位)：

$$L_{fs}=32.45+20{\log }_{10}(f)+20{\log }_{10}(d)$$

其中$f$为频率(MHz)，$d$为距离(km)。

2.2 大尺度衰落模型

2.2.1 对数距离路径损耗模型

实际环境中，路径损耗通常遵循对数距离规律：

$$PL(d)=PL(d_0)+10n{\log }_{10}\left(\frac{d}{d_0}\right)+X_{\sigma }$$

其中：

• $PL(d_0)$：参考距离$d_0$处的路径损耗(dB)
• $n$：路径损耗指数(自由空间$n=2$，城市环境$n=2.7\sim 4.0$)
• $X_{\sigma }$：阴影衰落，服从零均值高斯分布，标准差为$\sigma$

2.2.2 Okumura-Hata模型

Okumura-Hata模型是城市宏蜂窝最常用的经验模型之一，适用于150MHz-1500MHz频段：

城市地区：

$$L_{urban}=69.55+26.16{\log }_{10}(f)-13.82{\log }_{10}(h_b)-a(h_m)+[44.9-6.55{\log }_{10}(h_b)]{\log }_{10}(d)$$

其中：

• $f$：频率(MHz)
• $h_b$：基站天线高度(m)，30-200m
• $h_m$：移动台天线高度(m)，1-10m
• $d$：距离(km)，1-20km
• $a(h_m)$：移动台天线高度修正因子

对于中小城市：

$$a(h_m)=[1.1{\log }_{10}(f)-0.7]h_m-[1.56{\log }_{10}(f)-0.8]$$

对于大城市：

$$a(h_m)=8.29[{\log }_{10}(1.54h_m){]}^2-1.1,f\le 200\text{MHz}$$

$$a(h_m)=3.2[{\log }_{10}(11.75h_m){]}^2-4.97,f\ge 400\text{MHz}$$

郊区：

$$L_{suburban}=L_{urban}-2[{\log }_{10}(f/28){]}^2-5.4$$

开阔地：

$$L_{open}=L_{urban}-4.78[{\log }_{10}(f){]}^2+18.33{\log }_{10}(f)-40.94$$

2.2.3 COST-231 Hata模型

COST-231 Hata模型将频率范围扩展到2GHz：

$$L=46.3+33.9{\log }_{10}(f)-13.82{\log }_{10}(h_b)-a(h_m)+[44.9-6.55{\log }_{10}(h_b)]{\log }_{10}(d)+C_M$$

其中$C_M$为修正因子：

• 中等城市和郊区：$C_M=0$ dB
• 大都市中心：$C_M=3$ dB

2.3 小尺度衰落模型

2.3.1 瑞利衰落

当不存在直射路径(LOS)时，接收信号包络服从瑞利分布：

$$f(r)=\frac{r}{{\sigma }^2}\exp \left(-\frac{r^2}{2{\sigma }^2}\right),r\ge 0$$

其中${\sigma }^2$是多径分量的平均功率。

2.3.2 莱斯衰落

当存在强直射路径时，接收信号包络服从莱斯分布：

$$f(r)=\frac{r}{{\sigma }^2}\exp \left(-\frac{r^2+A^2}{2{\sigma }^2}\right)I_0\left(\frac{rA}{{\sigma }^2}\right),r\ge 0$$

其中：

• $A$：直射路径幅度
• $I_0$：第一类零阶修正贝塞尔函数
• $K=A^2/(2{\sigma }^2)$：莱斯K因子(dB)

2.3.3 多普勒频移

移动引起的多普勒频移：

$$f_d=\frac{v}{\lambda }\cos \theta =f_c\frac{v}{c}\cos \theta$$

其中：

• $v$：移动速度(m/s)
• $\theta$：入射波与移动方向的夹角
• $f_c$：载波频率(Hz)
• $c$：光速(m/s)

最大多普勒频移：

$$f_{d,max}=\frac{v}{\lambda }=f_c\frac{v}{c}$$

三、城市电波传播建模

3.1 城市传播环境特征

城市环境电波传播的主要特征包括：

1. 多径传播：信号经建筑物反射、散射后沿多条路径到达接收机
2. 阴影衰落：大型建筑物遮挡造成的信号衰减
3. 穿透损耗：信号穿透建筑物外墙、窗户等结构产生的损耗
4. 街道峡谷效应：高楼之间的街道形成波导-like传播

3.2 射线追踪方法

射线追踪(Ray Tracing)是基于几何光学(GO)和几何绕射理论(GTD/UTD)的确定性传播预测方法。

3.2.1 基本射线类型

直射射线(LOS)：

• 发射机与接收机之间无遮挡
• 场强计算：$E_{LOS}=E_0\frac{e^{-jkr}}{r}$

反射射线：

• 遵循反射定律：入射角等于反射角
• 反射系数：

对于垂直极化：

$$R_{\perp }=\frac{\cos {\theta }_i-\sqrt{{\epsilon }_r-{\sin }^2{\theta }_i}}{\cos {\theta }_i+\sqrt{{\epsilon }_r-{\sin }^2{\theta }_i}}$$

对于水平极化：

$$R_{\parallel }=\frac{{\epsilon }_r\cos {\theta }_i-\sqrt{{\epsilon }_r-{\sin }^2{\theta }_i}}{{\epsilon }_r\cos {\theta }_i+\sqrt{{\epsilon }_r-{\sin }^2{\theta }_i}}$$

其中${\epsilon }_r={\epsilon }_r^{\prime }-j\frac{\sigma }{\omega {\epsilon }_0}$为复介电常数。

绕射射线:

• 当射线遇到建筑物边缘或角落时发生绕射
• 使用UTD计算绕射系数
• 单刃绕射损耗近似：

$$L_d=-20{\log }_{10}\left|\frac{1}{2}-\frac{1}{2}\sqrt{\frac{1}{2}}\left(\frac{1}{2}+C(v)+S(v)\right)\right|$$

其中$C(v)$和$S(v)$为Fresnel积分，$v$为绕射参数：

$$v=h\sqrt{\frac{2(d_1+d_2)}{\lambda d_1{d}_2}}$$

散射射线:

• 粗糙表面或不规则结构引起的散射
• 散射场强通常远小于反射和绕射

3.2.2 射线追踪算法实现

射线追踪的基本步骤：

1. 场景建模：建立建筑物的三维几何模型
2. 射线发射：从发射点发射大量射线(通常采用镜像法或射线发射法)
3. 射线追踪：计算射线与场景中物体的交点
4. 场强计算：累加所有到达接收机的射线贡献
5. 路径损耗：转换为路径损耗或接收功率

加速技术：

• 空间划分(八叉树、KD树)
• 并行计算(GPU加速)
• 射线束追踪(Beam Tracing)

3.3 主导路径模型

主导路径模型(Dominant Path Model, DPM)是一种简化的射线追踪方法，只考虑对总场贡献最大的几条路径。

主导路径识别：

1. 计算直射路径(如果存在)
2. 识别主要的反射面(地面、大型建筑物墙面)
3. 计算关键绕射点(建筑物边缘、屋顶)
4. 评估每条路径的损耗，选择损耗最小的几条

路径损耗计算：

$$L_{total}=L_{fs}+\sum _i{L}_{reflection,i}+\sum _j{L}_{diffraction,j}+L_{penetration}$$

3.4 城市微蜂窝模型

3.4.1 街道峡谷模型

街道峡谷中的传播可以用波导模型近似：

$$PL=-10{\log }_{10}{\left(\frac{\lambda }{2\pi d}\right)}^2+10{\log }_{10}{\left(\frac{d}{d_0}\right)}^n+L_{excess}$$

其中$n$通常为4-6，$L_{excess}$为额外损耗。

COST-231 Walfisch-Ikegami模型：

适用于城市微蜂窝，频率800MHz-2GHz：

$$L=L_0+L_{rts}+L_{msd}$$

其中：

• $L_0$：自由空间损耗
• $L_{rts}$：屋顶到街道的绕射和散射损耗
• $L_{msd}$：多屏绕射损耗

屋顶到街道损耗：

$$L_{rts}=-16.9-10{\log }_{10}(w)+10{\log }_{10}(f)+20{\log }_{10}(\Delta h_m)+L_{ori}$$

其中：

• $w$：街道宽度(m)
• $\Delta h_m=h_{roof}-h_m$：建筑物高度与移动台高度差
• $L_{ori}$：方向损耗

多屏绕射损耗：

$$L_{msd}=L_{bsh}+k_a+k_d{\log }_{10}(d)+k_f{\log }_{10}(f)-9{\log }_{10}(b)$$

其中$b$为建筑物间距。

四、室内电波传播建模

4.1 室内传播环境特征

室内环境相比城市环境更加复杂：

1. 密集的多径：墙壁、地板、天花板、家具等造成大量反射
2. 穿透损耗变化大：不同材料(混凝土、砖、石膏板、玻璃)的损耗差异显著
3. 楼层间传播：信号通过楼层间的穿透或楼梯/电梯传播
4. 人员活动影响：人体遮挡和移动引起信号波动

4.2 室内路径损耗模型

4.2.1 对数距离模型(室内)

$$PL(d)=PL(d_0)+10n_{SF}{\log }_{10}\left(\frac{d}{d_0}\right)+FAF+X_{\sigma }$$

其中：

• $n_{SF}$：同楼层路径损耗指数(1.6-3.3)
• $FAF$：楼层衰减因子(Floor Attenuation Factor)
• $X_{\sigma }$：阴影衰落(3-6dB)

典型参数值：

环境类型	路径损耗指数$n$	标准差$\sigma$(dB)
自由空间	2	0
走廊	1.6-1.8	3-5
办公室(软隔断)	2.2-2.8	4-6
办公室(硬隔断)	2.8-3.3	5-7
工厂	2.0-3.0	5-8

4.2.2 楼层衰减因子(FAF)

$$FAF=\alpha \cdot N_F$$

其中$N_F$为楼层数，$\alpha$为每层衰减(通常10-20dB/层)。

4.2.3 穿透损耗

不同建筑材料的穿透损耗：

材料	1GHz损耗(dB)	2.4GHz损耗(dB)	5.8GHz损耗(dB)
混凝土墙(20cm)	10-15	15-20	20-25
砖墙	5-8	8-12	12-18
石膏板	2-4	3-5	5-8
玻璃窗	2-3	3-5	5-8
金属门	10-15	15-20	20-25
电梯	25-35	30-40	35-45

4.3 室内射线追踪

室内射线追踪需要考虑更多细节：

1. 精细几何建模：包括墙壁厚度、门窗位置、家具布局
2. 材料属性：不同表面的介电常数和电导率
3. 多次反射：通常考虑3-6次反射
4. 绕射：门缝、墙角等边缘绕射

加速技术：

• 可见性预处理
• 射线分类(直达、一次反射、二次反射等)
• 空间加速结构(BSP树、网格划分)

4.4 室内统计模型

4.4.1 Saleh-Valenzuela模型

室内多径信道的经典统计模型：

多径到达时间：

• 射线簇(Cluster)按泊松过程到达，速率$\Lambda$
• 簇内射线按泊松过程到达，速率$\lambda$

功率延迟分布：

$$P(\tau )=\sum _{l=0}^{\infty }\sum _{k=0}^{\infty }{P}_{k,l}\delta (\tau -T_l-{\tau }_{k,l})$$

其中：

• $T_l$：第$l$个簇的到达时间
• ${\tau }_{k,l}$：第$l$个簇中第$k$条射线的相对到达时间
• $P_{k,l}$：对应射线的功率

功率衰减：

• 簇内功率按指数衰减：$\overline{{\beta }^2}(T_l,{\tau }_{k,l})=\overline{{\beta }^2(0,0)}{e}^{-T_l/\Gamma }{e}^{-{\tau }_{k,l}/\gamma }$
• $\Gamma$：簇衰减时间常数
• $\gamma$：射线衰减时间常数

4.4.2 IEEE 802.11信道模型

IEEE 802.11标准定义了多种室内信道模型：

Model A(办公室环境，LOS)：

• RMS时延扩展：50ns
• 最大时延：100ns

Model B(办公室环境，NLOS)：

• RMS时延扩展：100ns
• 最大时延：200ns

Model C(大型开放空间/室内NLOS)：

• RMS时延扩展：150ns
• 最大时延：500ns

Model D(大型开放空间，NLOS)：

• RMS时延扩展：140ns
• 最大时延：390ns

Model E(大型开放空间，NLOS)：

• RMS时延扩展：250ns
• 最大时延：730ns

五、5G/6G毫米波传播特性

5.1 毫米波传播特点

毫米波(24GHz以上)传播与低频段有显著差异：

1. 高路径损耗：自由空间损耗随频率平方增加
2. 强方向性：需要使用高增益定向天线
3. 大气吸收：氧气(60GHz)和水蒸气(183GHz)吸收峰
4. 穿透损耗大：难以穿透建筑物和人体
5. 雨衰严重：高频段受降雨影响显著

5.2 毫米波信道模型

5.2.1 3GPP TR 38.901模型

3GPP为5G NR定义的信道模型，支持0.5-100GHz：

路径损耗模型：

UMi(城市微蜂窝)：

$$P{L}_{UMi}=32.4+21{\log }_{10}(d_{3D})+20{\log }_{10}(f_c),\text{LOS}$$

$$P{L}_{UMi}=35.3{\log }_{10}(d_{3D})+22.4+21.3{\log }_{10}(f_c)-0.3(h_{UT}-1.5),\text{NLOS}$$

UMa(城市宏蜂窝)：

$$P{L}_{UMa}=28.0+22{\log }_{10}(d_{3D})+20{\log }_{10}(f_c),\text{LOS}$$

$$P{L}_{UMa}=13.54+39.08{\log }_{10}(d_{3D})+20{\log }_{10}(f_c)-0.6(h_{UT}-1.5),\text{NLOS}$$

RMa(农村宏蜂窝)：

$$P{L}_{RMa}=31.4+21.5{\log }_{10}(d_{3D})+19.1{\log }_{10}(f_c),\text{LOS}$$

$$P{L}_{RMa}=max(P{L}_{RMa-LOS},33.83+31.94{\log }_{10}(d_{3D})+20{\log }_{10}(f_c)),\text{NLOS}$$

InH(室内热点)：

$$P{L}_{InH}=32.4+17.3{\log }_{10}(d_{3D})+20{\log }_{10}(f_c),\text{LOS}$$

$$P{L}_{InH}=33.6+21.9{\log }_{10}(d_{3D})+20{\log }_{10}(f_c),\text{NLOS}$$

其中：

• $d_{3D}$：三维距离(m)
• $f_c$：载波频率(GHz)
• $h_{UT}$：用户终端高度(m)

5.2.2 毫米波穿透损耗

室外到室内(O2I)：

$$P{L}_{O2I}=P{L}_{outdoor}+P{L}_{tw}+P{L}_{in}$$

其中：

• $P{L}_{tw}$：穿透外墙损耗(5-20dB，取决于材料和频率)
• $P{L}_{in}$：室内传播损耗

标准穿透损耗值(3GPP)：

场景	低频(<6GHz)	高频(>6GHz)
标准玻璃	3.5 dB	6.5 dB
低辐射玻璃	15 dB	30 dB
混凝土墙	15 dB	30 dB
木门	4 dB	8 dB

5.3 波束管理与覆盖

毫米波系统依赖波束赋形来克服高路径损耗：

波束扫描：

• 发射端和接收端通过扫描窄波束寻找最佳路径
• 波束宽度通常为10°-30°

波束跟踪：

• 用户移动时动态调整波束方向
• 需要快速跟踪算法

波束切换：

• 当主波束被遮挡时切换到备用波束
• 切换时间要求<10ms

六、Python仿真实现

6.1 仿真1：自由空间与基本路径损耗模型

"""
仿真1：自由空间与基本路径损耗模型
比较不同路径损耗模型的预测结果
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import os

# 设置Agg后端，不弹出窗口
import matplotlib
matplotlib.use('Agg')

# 创建输出目录
output_dir = os.path.dirname(os.path.abspath(__file__))

print("=" * 60)
print("仿真1: 自由空间与基本路径损耗模型")
print("=" * 60)

# 参数设置
frequencies = [900e6, 1800e6, 2100e6, 2600e6, 3500e6, 28e9]  # Hz
freq_labels = ['900MHz', '1.8GHz', '2.1GHz', '2.6GHz', '3.5GHz', '28GHz']
distances = np.linspace(1, 1000, 1000)  # m

# 自由空间路径损耗
def free_space_path_loss(f, d):
    """计算自由空间路径损耗 (dB)"""
    c = 3e8  # 光速
    lambda_ = c / f
    # 确保距离大于0
    d = np.maximum(d, 0.001)
    fspl = 20 * np.log10(d) + 20 * np.log10(f) + 20 * np.log10(4*np.pi/c)
    return fspl

# 对数距离路径损耗模型
def log_distance_path_loss(f, d, d0, n, sigma=0):
    """
    对数距离路径损耗模型
    f: 频率(Hz)
    d: 距离(m)
    d0: 参考距离(m)
    n: 路径损耗指数
    sigma: 阴影衰落标准差(dB)
    """
    d = np.maximum(d, 0.001)
    pl0 = free_space_path_loss(f, d0)
    pl = pl0 + 10 * n * np.log10(d / d0)
    if sigma > 0:
        pl += np.random.normal(0, sigma, size=np.shape(d))
    return pl

# Okumura-Hata模型 (城市)
def okumura_hata_urban(f, d, hb, hm):
    """
    Okumura-Hata城市模型
    f: 频率(MHz)
    d: 距离(km)
    hb: 基站天线高度(m)
    hm: 移动台天线高度(m)
    """
    # 频率范围检查
    if f < 150 or f > 1500:
        print(f"警告: 频率{f}MHz超出Okumura-Hata模型适用范围(150-1500MHz)")
    
    # 移动台天线高度修正因子 (中小城市)
    a_hm = (1.1 * np.log10(f) - 0.7) * hm - (1.56 * np.log10(f) - 0.8)
    
    # 路径损耗
    L = 69.55 + 26.16 * np.log10(f) - 13.82 * np.log10(hb) - a_hm
    L += (44.9 - 6.55 * np.log10(hb)) * np.log10(d)
    
    return L

# COST-231 Hata模型
def cost231_hata(f, d, hb, hm, urban_type='medium'):
    """
    COST-231 Hata模型
    f: 频率(MHz)
    d: 距离(km)
    hb: 基站天线高度(m)
    hm: 移动台天线高度(m)
    urban_type: 'medium'或'metropolitan'
    """
    # 频率范围检查
    if f < 1500 or f > 2000:
        print(f"警告: 频率{f}MHz超出COST-231模型最佳适用范围(1500-2000MHz)")
    
    # 移动台天线高度修正因子
    a_hm = (1.1 * np.log10(f) - 0.7) * hm - (1.56 * np.log10(f) - 0.8)
    
    # 大都市修正因子
    C_M = 3 if urban_type == 'metropolitan' else 0
    
    # 路径损耗
    L = 46.3 + 33.9 * np.log10(f) - 13.82 * np.log10(hb) - a_hm
    L += (44.9 - 6.55 * np.log10(hb)) * np.log10(d) + C_M
    
    return L

# 计算各种模型的路径损耗
print("\n计算路径损耗模型...")

# 选择代表性频率进行分析
f_analysis = 1800e6  # 1.8GHz
f_mhz = f_analysis / 1e6
d_km = distances / 1000

# 1. 自由空间路径损耗
fspl = free_space_path_loss(f_analysis, distances)

# 2. 对数距离模型 (不同路径损耗指数)
pl_n2 = log_distance_path_loss(f_analysis, distances, d0=100, n=2.0)
pl_n3 = log_distance_path_loss(f_analysis, distances, d0=100, n=3.0)
pl_n4 = log_distance_path_loss(f_analysis, distances, d0=100, n=4.0)

# 3. Okumura-Hata模型
hata_urban = okumura_hata_urban(f_mhz, d_km, hb=30, hm=1.5)

# 4. COST-231 Hata模型
cost231_medium = cost231_hata(f_mhz, d_km, hb=30, hm=1.5, urban_type='medium')
cost231_metro = cost231_hata(f_mhz, d_km, hb=30, hm=1.5, urban_type='metropolitan')

print("  完成路径损耗计算")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 不同频率的自由空间路径损耗
ax = axes[0, 0]
colors = plt.cm.viridis(np.linspace(0, 1, len(frequencies)))
for f, label, color in zip(frequencies, freq_labels, colors):
    pl = free_space_path_loss(f, distances)
    ax.plot(distances/1000, pl, color=color, linewidth=2, label=label)
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Path Loss (dB)')
ax.set_title('Free Space Path Loss vs Frequency')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 1)

# 图2: 不同路径损耗模型的比较
ax = axes[0, 1]
ax.plot(distances/1000, fspl, 'k-', linewidth=2, label='Free Space (n=2)')
ax.plot(distances/1000, pl_n3, 'b--', linewidth=2, label='Log-Distance (n=3)')
ax.plot(distances/1000, pl_n4, 'r-.', linewidth=2, label='Log-Distance (n=4)')
ax.plot(distances/1000, hata_urban, 'g:', linewidth=2, label='Okumura-Hata')
ax.plot(distances/1000, cost231_medium, 'm-', linewidth=2, label='COST-231')
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Path Loss (dB)')
ax.set_title(f'Path Loss Models Comparison @ {f_mhz:.0f}MHz')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 1)

# 图3: 阴影衰落的影响
ax = axes[1, 0]
np.random.seed(42)
pl_shadow = log_distance_path_loss(f_analysis, distances, d0=100, n=3.0, sigma=6)
ax.plot(distances/1000, pl_n3, 'b-', linewidth=2, label='Median (σ=0)')
ax.fill_between(distances/1000, pl_n3 - 6, pl_n3 + 6, alpha=0.3, label='±1σ')
ax.fill_between(distances/1000, pl_n3 - 12, pl_n3 + 12, alpha=0.2, label='±2σ')
ax.plot(distances/1000, pl_shadow, 'r-', alpha=0.5, linewidth=1, label='Sample with σ=6dB')
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Path Loss (dB)')
ax.set_title('Shadow Fading Effect (Log-Normal)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 1)

# 图4: 频率对路径损耗的影响
ax = axes[1, 1]
d_fixed = 500  # m
f_range = np.linspace(100e6, 6e9, 100)
fspl_freq = free_space_path_loss(f_range, d_fixed)
# 转换为MHz用于显示
f_mhz_range = f_range / 1e6
ax.plot(f_mhz_range, fspl_freq, 'b-', linewidth=2)
ax.axvline(x=700, color='r', linestyle='--', alpha=0.5, label='LTE Band')
ax.axvline(x=3500, color='g', linestyle='--', alpha=0.5, label='5G Sub-6')
ax.axvline(x=28000, color='m', linestyle='--', alpha=0.5, label='5G mmWave')
ax.set_xlabel('Frequency (MHz)')
ax.set_ylabel('Path Loss (dB)')
ax.set_title(f'Path Loss vs Frequency @ d={d_fixed}m')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
output_path = os.path.join(output_dir, 'path_loss_models.png')
plt.savefig(output_path, dpi=150, bbox_inches='tight')
plt.close()
print(f"Saved: {output_path}")

# 打印关键数值
print("\n关键数值结果:")
print(f"  1.8GHz, 1km处:")
print(f"    自由空间损耗: {free_space_path_loss(1800e6, 1000):.1f} dB")
print(f"    Okumura-Hata(城市): {okumura_hata_urban(1800, 1, 30, 1.5):.1f} dB")
print(f"    COST-231(中等城市): {cost231_hata(1800, 1, 30, 1.5, 'medium'):.1f} dB")
print(f"    COST-231(大都市): {cost231_hata(1800, 1, 30, 1.5, 'metropolitan'):.1f} dB")

6.2 仿真2：城市街道峡谷传播模型

"""
仿真2：城市街道峡谷传播模型
模拟城市街道峡谷中的电波传播特性
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch
import os

matplotlib.use('Agg')

print("\n" + "=" * 60)
print("仿真2: 城市街道峡谷传播模型")
print("=" * 60)

# COST-231 Walfisch-Ikegami模型
def walfisch_ikegami(f, d, hb, hm, hroof, w, phi, environment='urban'):
    """
    COST-231 Walfisch-Ikegami模型
    f: 频率(MHz)
    d: 距离(km)
    hb: 基站天线高度(m)
    hm: 移动台天线高度(m)
    hroof: 建筑物平均高度(m)
    w: 街道宽度(m)
    phi: 街道与发射方向的夹角(度)
    environment: 'urban', 'suburban', 'rural'
    """
    # 自由空间损耗
    L0 = 32.45 + 20*np.log10(f) + 20*np.log10(d*1000)  # d转换为m
    
    # 检查是否为LOS条件
    delta_hb = hb - hroof
    delta_hm = hroof - hm
    
    # LOS条件简化处理：如果基站高于屋顶且距离较近
    if delta_hb > 0 and d < 0.5:
        return L0  # 近似LOS
    
    # 屋顶到街道损耗
    if phi < 35:
        Lori = -10 + 0.354*phi
    elif phi < 55:
        Lori = 2.5 + 0.075*(phi - 35)
    else:
        Lori = 4.0 - 0.114*(phi - 55)
    
    Lrts = -8.2 - 10*np.log10(w) + 10*np.log10(f) + 20*np.log10(delta_hm) + Lori
    
    # 多屏绕射损耗
    if hb > hroof:
        Lbsh = -18*np.log10(1 + delta_hb)
        ka = 54
        kd = 18
    else:
        Lbsh = 0
        ka = 54 - 0.8*delta_hb
        kd = 18 - 15*delta_hb/hroof
    
    # 频率修正
    kf = -4 + 0.7*(f/925 - 1)
    
    # 建筑物间距(假设)
    b = 50  # m
    
    Lmsd = Lbsh + ka + kd*np.log10(d) + kf*np.log10(f) - 9*np.log10(b)
    
    # 总损耗
    L = L0 + Lrts + Lmsd
    
    # 郊区修正
    if environment == 'suburban':
        L -= 2*(np.log10(f/28))**2 - 5.4
    elif environment == 'rural':
        L -= 4.78*(np.log10(f))**2 - 18.33*np.log10(f) - 40.94
    
    return L

# 街道峡谷场景参数
f = 1800  # MHz
d_range = np.linspace(0.1, 2, 100)  # km
hb = 30  # m (基站高度)
hm = 1.5  # m (移动台高度)
hroof = 20  # m (建筑物高度)
w = 20  # m (街道宽度)

# 不同街道角度
phi_values = [0, 45, 90]  # 度
phi_labels = ['0° (平行)', '45°', '90° (垂直)']

print("\n计算街道峡谷传播损耗...")

# 计算不同角度的路径损耗
losses = {}
for phi, label in zip(phi_values, phi_labels):
    losses[label] = [walfisch_ikegami(f, d, hb, hm, hroof, w, phi) for d in d_range]

print("  完成计算")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 不同街道角度的路径损耗
ax = axes[0, 0]
colors = ['b', 'g', 'r']
for (label, loss), color in zip(losses.items(), colors):
    ax.plot(d_range, loss, color=color, linewidth=2, label=label)
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Path Loss (dB)')
ax.set_title(f'Street Canyon Propagation @ {f}MHz\n(hb={hb}m, hroof={hroof}m, w={w}m)')
ax.legend()
ax.grid(True, alpha=0.3)

# 图2: 建筑物高度的影响
ax = axes[0, 1]
hroof_values = [10, 15, 20, 25, 30]
colors = plt.cm.viridis(np.linspace(0, 1, len(hroof_values)))
for hroof_test, color in zip(hroof_values, colors):
    loss_hroof = [walfisch_ikegami(f, d, hb, hm, hroof_test, w, 90) for d in d_range]
    ax.plot(d_range, loss_hroof, color=color, linewidth=2, label=f'hroof={hroof_test}m')
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Path Loss (dB)')
ax.set_title('Effect of Building Height')
ax.legend()
ax.grid(True, alpha=0.3)

# 图3: 街道宽度的影响
ax = axes[1, 0]
w_values = [10, 20, 30, 40, 50]
colors = plt.cm.plasma(np.linspace(0, 1, len(w_values)))
for w_test, color in zip(w_values, colors):
    loss_w = [walfisch_ikegami(f, d, hb, hm, hroof, w_test, 90) for d in d_range]
    ax.plot(d_range, loss_w, color=color, linewidth=2, label=f'w={w_test}m')
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Path Loss (dB)')
ax.set_title('Effect of Street Width')
ax.legend()
ax.grid(True, alpha=0.3)

# 图4: 场景示意图
ax = axes[1, 1]
ax.set_xlim(-50, 150)
ax.set_ylim(-5, 40)

# 绘制建筑物
building1 = Rectangle((-40, 0), 40, hroof, facecolor='gray', edgecolor='black')
building2 = Rectangle((w, 0), 40, hroof, facecolor='gray', edgecolor='black')
ax.add_patch(building1)
ax.add_patch(building2)

# 绘制基站
ax.plot([0], [hb], 'ro', markersize=12, label='Base Station')
ax.plot([0, 0], [0, hb], 'r--', alpha=0.5)

# 绘制移动台
ax.plot([100], [hm], 'bs', markersize=10, label='Mobile')
ax.plot([100, 100], [0, hm], 'b--', alpha=0.5)

# 绘制射线
ax.plot([0, 100], [hb, hm], 'g-', linewidth=2, alpha=0.6, label='Direct path')
ax.plot([0, w/2, 100], [hb, hroof, hm], 'm:', linewidth=2, alpha=0.6, label='Diffracted path')

ax.set_xlabel('Distance (m)')
ax.set_ylabel('Height (m)')
ax.set_title('Street Canyon Scenario')
ax.legend(loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_aspect('equal', adjustable='box')

plt.tight_layout()
output_path = os.path.join(output_dir, 'street_canyon_propagation.png')
plt.savefig(output_path, dpi=150, bbox_inches='tight')
plt.close()
print(f"Saved: {output_path}")

print("\n关键结果:")
print(f"  1.8GHz, 1km, 垂直街道:")
print(f"    路径损耗: {walfisch_ikegami(1800, 1, 30, 1.5, 20, 20, 90):.1f} dB")

6.3 仿真3：室内传播模型与穿透损耗

"""
仿真3：室内传播模型与穿透损耗
分析室内环境中的电波传播特性
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import os

matplotlib.use('Agg')

print("\n" + "=" * 60)
print("仿真3: 室内传播模型与穿透损耗")
print("=" * 60)

# 室内路径损耗模型
def indoor_path_loss(f, d, n, d0=1, sigma=0, FAF=0):
    """
    室内对数距离路径损耗模型
    f: 频率(GHz)
    d: 距离(m)
    n: 路径损耗指数
    d0: 参考距离(m)
    sigma: 阴影衰落标准差(dB)
    FAF: 楼层衰减因子(dB)
    """
    # 参考距离处的自由空间损耗
    c = 3e8
    lambda_ = c / (f * 1e9)
    PL0 = 20 * np.log10(4 * np.pi * d0 / lambda_)
    
    # 对数距离损耗
    d = np.maximum(d, 0.1)
    PL = PL0 + 10 * n * np.log10(d / d0) + FAF
    
    if sigma > 0:
        PL += np.random.normal(0, sigma)
    
    return PL

# 穿透损耗数据 (dB)
penetration_loss = {
    'Concrete (20cm)': {'1GHz': 12, '2.4GHz': 18, '5.8GHz': 25},
    'Brick Wall': {'1GHz': 6, '2.4GHz': 10, '5.8GHz': 15},
    'Drywall': {'1GHz': 3, '2.4GHz': 5, '5.8GHz': 8},
    'Glass Window': {'1GHz': 2, '2.4GHz': 4, '5.8GHz': 7},
    'Wooden Door': {'1GHz': 3, '2.4GHz': 5, '5.8GHz': 8},
    'Metal Door': {'1GHz': 15, '2.4GHz': 25, '5.8GHz': 35},
}

# 不同环境的路径损耗指数
environment_params = {
    'Free Space': {'n': 2.0, 'sigma': 0},
    'Corridor': {'n': 1.7, 'sigma': 3.0},
    'Office (Soft Partition)': {'n': 2.4, 'sigma': 5.0},
    'Office (Hard Partition)': {'n': 3.0, 'sigma': 6.0},
    'Factory': {'n': 2.5, 'sigma': 7.0},
}

print("\n计算室内传播损耗...")

# 频率和距离范围
frequencies = [0.9, 2.4, 5.8]  # GHz
distances = np.linspace(1, 100, 100)  # m

# 计算不同环境的损耗
losses_env = {}
for env, params in environment_params.items():
    losses_env[env] = [indoor_path_loss(2.4, d, params['n'], sigma=params['sigma']) 
                       for d in distances]

print("  完成计算")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 不同室内环境的路径损耗
ax = axes[0, 0]
colors = plt.cm.tab10(np.linspace(0, 1, len(environment_params)))
for (env, loss), color in zip(losses_env.items(), colors):
    ax.plot(distances, loss, color=color, linewidth=2, label=env)
ax.set_xlabel('Distance (m)')
ax.set_ylabel('Path Loss (dB)')
ax.set_title('Indoor Path Loss by Environment @ 2.4GHz')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 100)

# 图2: 不同频率的室内损耗
ax = axes[0, 1]
np.random.seed(42)
colors = ['b', 'g', 'r']
for f, color in zip(frequencies, colors):
    # 办公室环境(软隔断)
    loss_f = [indoor_path_loss(f, d, 2.4, sigma=5) for d in distances]
    ax.plot(distances, loss_f, color=color, linewidth=2, label=f'{f}GHz')
ax.set_xlabel('Distance (m)')
ax.set_ylabel('Path Loss (dB)')
ax.set_title('Indoor Path Loss vs Frequency\n(Office, Soft Partition)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 100)

# 图3: 穿透损耗比较
ax = axes[1, 0]
materials = list(penetration_loss.keys())
freqs_plot = ['1GHz', '2.4GHz', '5.8GHz']
x = np.arange(len(materials))
width = 0.25

for i, freq in enumerate(freqs_plot):
    losses_material = [penetration_loss[m][freq] for m in materials]
    ax.bar(x + i*width, losses_material, width, label=freq)

ax.set_xlabel('Material')
ax.set_ylabel('Penetration Loss (dB)')
ax.set_title('Penetration Loss by Material and Frequency')
ax.set_xticks(x + width)
ax.set_xticklabels(materials, rotation=45, ha='right', fontsize=8)
ax.legend()
ax.grid(True, alpha=0.3, axis='y')

# 图4: 楼层衰减因子
ax = axes[1, 1]
floors = np.arange(0, 6)
FAF_values = [0, 15, 28, 38, 45, 50]  # dB (典型值)
ax.bar(floors, FAF_values, color='steelblue', alpha=0.7)
ax.set_xlabel('Number of Floors')
ax.set_ylabel('Floor Attenuation Factor (dB)')
ax.set_title('Floor Attenuation Factor (FAF)\n@ 2.4GHz')
ax.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
output_path = os.path.join(output_dir, 'indoor_propagation.png')
plt.savefig(output_path, dpi=150, bbox_inches='tight')
plt.close()
print(f"Saved: {output_path}")

print("\n关键结果:")
print("  2.4GHz, 50m处:")
for env, params in environment_params.items():
    loss = indoor_path_loss(2.4, 50, params['n'])
    print(f"    {env}: {loss:.1f} dB")

6.4 仿真4：射线追踪基础实现

"""
仿真4：射线追踪基础实现
实现简单的2D射线追踪算法
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, Circle, FancyArrowPatch
import os

matplotlib.use('Agg')

print("\n" + "=" * 60)
print("仿真4: 射线追踪基础实现")
print("=" * 60)

# 简单的射线类
class Ray:
    def __init__(self, origin, direction):
        self.origin = np.array(origin)
        self.direction = np.array(direction) / np.linalg.norm(direction)
        self.path = [self.origin.copy()]
        self.reflections = 0
    
    def point_at(self, t):
        return self.origin + t * self.direction
    
    def reflect(self, normal):
        """计算反射方向"""
        n = np.array(normal) / np.linalg.norm(normal)
        self.direction = self.direction - 2 * np.dot(self.direction, n) * n
        self.reflections += 1

# 墙壁类
class Wall:
    def __init__(self, x1, y1, x2, y2, material='concrete'):
        self.p1 = np.array([x1, y1])
        self.p2 = np.array([x2, y2])
        self.material = material
        # 计算法向量(指向外侧)
        dx = x2 - x1
        dy = y2 - y1
        self.normal = np.array([-dy, dx])
        self.normal = self.normal / np.linalg.norm(self.normal)
    
    def intersect(self, ray):
        """计算射线与线段的交点"""
        # 线段参数方程: P = p1 + t*(p2-p1)
        # 射线参数方程: P = origin + s*direction
        v1 = ray.origin - self.p1
        v2 = self.p2 - self.p1
        v3 = np.array([-ray.direction[1], ray.direction[0]])
        
        denom = np.dot(v2, v3)
        if abs(denom) < 1e-6:
            return None  # 平行
        
        t1 = np.cross(v2, v1) / denom
        t2 = np.dot(v1, v3) / denom
        
        if t1 >= 0 and 0 <= t2 <= 1:
            return ray.point_at(t1)
        return None

# 场景设置
print("\n设置射线追踪场景...")

# 创建简单的室内场景(房间)
walls = [
    Wall(0, 0, 20, 0, 'concrete'),      # 底墙
    Wall(20, 0, 20, 15, 'concrete'),    # 右墙
    Wall(20, 15, 0, 15, 'concrete'),    # 顶墙
    Wall(0, 15, 0, 0, 'concrete'),      # 左墙
    Wall(8, 0, 8, 6, 'drywall'),        # 内墙1
    Wall(12, 9, 12, 15, 'drywall'),     # 内墙2
]

# 发射机和接收机位置
tx_pos = np.array([5, 7.5])
rx_pos = np.array([15, 7.5])

# 反射系数(简化)
reflection_coeff = {
    'concrete': 0.7,
    'drywall': 0.5,
    'glass': 0.3,
}

print("  场景设置完成")
print(f"  发射机位置: ({tx_pos[0]}, {tx_pos[1]})")
print(f"  接收机位置: ({rx_pos[0]}, {rx_pos[1]})")

# 发射射线
print("\n发射射线并追踪...")
n_rays = 36  # 发射36条射线，间隔10度
rays = []
received_power = []

for i in range(n_rays):
    angle = 2 * np.pi * i / n_rays
    direction = np.array([np.cos(angle), np.sin(angle)])
    ray = Ray(tx_pos, direction)
    
    # 追踪射线(最多3次反射)
    max_reflections = 3
    current_pos = tx_pos.copy()
    
    for _ in range(max_reflections + 1):
        # 找到最近的交点
        closest_intersection = None
        closest_wall = None
        min_distance = float('inf')
        
        for wall in walls:
            intersection = wall.intersect(ray)
            if intersection is not None:
                dist = np.linalg.norm(intersection - current_pos)
                if dist > 0.01 and dist < min_distance:  # 避免重复检测同一点
                    min_distance = dist
                    closest_intersection = intersection
                    closest_wall = wall
        
        if closest_intersection is not None:
            # 记录路径点
            ray.path.append(closest_intersection.copy())
            
            # 检查是否到达接收机附近
            if np.linalg.norm(closest_intersection - rx_pos) < 0.5:
                received_power.append({
                    'ray_id': i,
                    'angle': np.rad2deg(angle),
                    'reflections': ray.reflections,
                    'path_length': sum([np.linalg.norm(ray.path[j+1] - ray.path[j]) 
                                       for j in range(len(ray.path)-1)]),
                    'material': closest_wall.material if closest_wall else 'none'
                })
                break
            
            # 反射
            if ray.reflections < max_reflections:
                ray.origin = closest_intersection
                # 确保法向量指向正确的方向
                if np.dot(ray.direction, closest_wall.normal) > 0:
                    wall_normal = -closest_wall.normal
                else:
                    wall_normal = closest_wall.normal
                ray.reflect(wall_normal)
                current_pos = closest_intersection.copy()
            else:
                break
        else:
            # 射线射出场景
            ray.path.append(ray.point_at(30))  # 延伸到远处
            break
    
    rays.append(ray)

print(f"  发射了 {n_rays} 条射线")
print(f"  到达接收机的射线: {len(received_power)} 条")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 图1: 射线追踪场景
ax = axes[0]

# 绘制墙壁
for wall in walls:
    color = 'gray' if wall.material == 'concrete' else 'lightblue'
    ax.plot([wall.p1[0], wall.p2[0]], [wall.p1[1], wall.p2[1]], 
            'k-', linewidth=3)
    ax.fill_between([wall.p1[0], wall.p2[0]], 
                    [wall.p1[1], wall.p2[1]], 
                    [wall.p1[1], wall.p2[1]], 
                    alpha=0.3, color=color)

# 绘制射线(只显示部分以避免混乱)
for i, ray in enumerate(rays):
    if len(ray.path) > 1:
        path_array = np.array(ray.path)
        alpha = 0.3 if len(received_power) == 0 or i not in [p['ray_id'] for p in received_power] else 0.8
        linewidth = 1 if alpha == 0.3 else 2
        color = 'blue' if alpha == 0.3 else 'red'
        ax.plot(path_array[:, 0], path_array[:, 1], 
                color=color, alpha=alpha, linewidth=linewidth)

# 绘制发射机和接收机
ax.plot(tx_pos[0], tx_pos[1], 'g^', markersize=15, label='Transmitter')
ax.plot(rx_pos[0], rx_pos[1], 'ro', markersize=12, label='Receiver')

ax.set_xlabel('X (m)')
ax.set_ylabel('Y (m)')
ax.set_title('2D Ray Tracing Simulation')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
ax.set_xlim(-1, 21)
ax.set_ylim(-1, 16)

# 图2: 到达接收机的射线分析
ax = axes[1]

if len(received_power) > 0:
    angles = [p['angle'] for p in received_power]
    path_lengths = [p['path_length'] for p in received_power]
    reflections = [p['reflections'] for p in received_power]
    
    # 极坐标图显示到达角度
    ax2 = plt.subplot(122, projection='polar')
    theta = np.deg2rad(angles)
    r = [20 - pl/2 for pl in path_lengths]  # 路径越短，距离中心越近
    colors = plt.cm.viridis(np.array(reflections) / max(reflections))
    ax2.scatter(theta, r, c=reflections, s=100, cmap='viridis', alpha=0.7)
    ax2.set_title('Arrival Angles at Receiver\n(color = # reflections)')
    ax2.set_ylim(0, 20)
else:
    ax.text(0.5, 0.5, 'No rays reached receiver\nin this configuration', 
            ha='center', va='center', transform=ax.transAxes, fontsize=14)
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)
    ax.axis('off')

plt.tight_layout()
output_path = os.path.join(output_dir, 'ray_tracing_basic.png')
plt.savefig(output_path, dpi=150, bbox_inches='tight')
plt.close()
print(f"Saved: {output_path}")

if len(received_power) > 0:
    print("\n到达接收机的射线信息:")
    for i, p in enumerate(received_power[:5]):  # 显示前5条
        print(f"  射线 {p['ray_id']}: 角度={p['angle']:.1f}°, "
              f"反射={p['reflections']}次, 路径长={p['path_length']:.2f}m")

6.5 仿真5：5G毫米波传播特性

"""
仿真5：5G毫米波传播特性
分析毫米波频段的传播特性，包括路径损耗、穿透损耗、雨衰等
"""

import numpy as np
import matplotlib.pyplot as plt
import os

matplotlib.use('Agg')

print("\n" + "=" * 60)
print("仿真5: 5G毫米波传播特性")
print("=" * 60)

# 3GPP TR 38.901路径损耗模型
def path_loss_3gpp_umi_los(f, d):
    """3GPP UMi街道峡谷LOS模型"""
    # f in GHz, d in m
    PL = 32.4 + 21*np.log10(d) + 20*np.log10(f)
    return PL

def path_loss_3gpp_umi_nlos(f, d, h_ut=1.5):
    """3GPP UMi街道峡谷NLOS模型"""
    PL = 35.3*np.log10(d) + 22.4 + 21.3*np.log10(f) - 0.3*(h_ut - 1.5)
    return PL

def path_loss_3gpp_uma_los(f, d):
    """3GPP UMa LOS模型"""
    PL = 28.0 + 22*np.log10(d) + 20*np.log10(f)
    return PL

def path_loss_3gpp_uma_nlos(f, d, h_ut=1.5):
    """3GPP UMa NLOS模型"""
    PL = 13.54 + 39.08*np.log10(d) + 20*np.log10(f) - 0.6*(h_ut - 1.5)
    return PL

def path_loss_3gpp_inh_los(f, d):
    """3GPP InH LOS模型"""
    PL = 32.4 + 17.3*np.log10(d) + 20*np.log10(f)
    return PL

def path_loss_3gpp_inh_nlos(f, d):
    """3GPP InH NLOS模型"""
    PL = 33.6 + 21.9*np.log10(d) + 20*np.log10(f)
    return PL

# 大气吸收损耗
def atmospheric_absorption(f, d):
    """
    计算大气吸收损耗
    f: 频率(GHz)
    d: 距离(km)
    返回: 损耗(dB)
    """
    # 简化的氧气和水蒸气吸收模型
    # 氧气吸收峰在60GHz
    # 水蒸气吸收峰在183GHz
    
    if 50 <= f <= 70:
        # 氧气吸收带
        alpha = 0.1 * np.exp(-((f - 60)**2) / 50)  # dB/km
    elif 170 <= f <= 200:
        # 水蒸气吸收带
        alpha = 0.5 * np.exp(-((f - 183)**2) / 200)  # dB/km
    else:
        # 其他频率
        alpha = 0.01 + 0.001 * f  # dB/km, 简化模型
    
    return alpha * d

# 雨衰模型 (简化ITU-R模型)
def rain_attenuation(f, d, rain_rate):
    """
    计算雨衰
    f: 频率(GHz)
    d: 距离(km)
    rain_rate: 降雨率(mm/h)
    返回: 损耗(dB)
"""
    # 简化模型
    if rain_rate <= 0:
        return 0
    
    # 频率相关性
    k = 0.0001 * f**1.5  # 简化系数
    alpha = 0.8 + 0.02 * f  # 简化指数
    
    # 特定衰减 (dB/km)
    gamma = k * rain_rate**alpha
    
    # 有效路径长度修正
    deff = d * (1 - np.exp(-d/5)) / (d/5) if d > 0 else 0
    
    return gamma * deff

# 穿透损耗 (3GPP标准值)
penetration_loss_3gpp = {
    'Standard Glass': {'LF': 3.5, 'HF': 6.5},      # dB
    'Low-E Glass': {'LF': 15, 'HF': 30},
    'Concrete Wall': {'LF': 15, 'HF': 30},
    'Wooden Door': {'LF': 4, 'HF': 8},
}

print("\n计算5G毫米波传播特性...")

# 频率范围
frequencies = [0.7, 1.8, 2.6, 3.5, 28, 39, 60]  # GHz
freq_labels = ['0.7GHz', '1.8GHz', '2.6GHz', '3.5GHz', '28GHz', '39GHz', '60GHz']

# 距离范围
distances = np.linspace(10, 500, 100)  # m

# 计算不同场景的路径损耗
print("  计算UMi场景...")
loss_umi_los = [path_loss_3gpp_umi_los(28, d) for d in distances]
loss_umi_nlos = [path_loss_3gpp_umi_nlos(28, d) for d in distances]

print("  计算UMa场景...")
loss_uma_los = [path_loss_3gpp_uma_los(28, d) for d in distances]
loss_uma_nlos = [path_loss_3gpp_uma_nlos(28, d) for d in distances]

print("  计算InH场景...")
loss_inh_los = [path_loss_3gpp_inh_los(28, d) for d in distances]
loss_inh_nlos = [path_loss_3gpp_inh_nlos(28, d) for d in distances]

print("  完成计算")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 不同场景的路径损耗比较
ax = axes[0, 0]
ax.plot(distances, loss_umi_los, 'b-', linewidth=2, label='UMi LOS')
ax.plot(distances, loss_umi_nlos, 'b--', linewidth=2, label='UMi NLOS')
ax.plot(distances, loss_uma_los, 'r-', linewidth=2, label='UMa LOS')
ax.plot(distances, loss_uma_nlos, 'r--', linewidth=2, label='UMa NLOS')
ax.plot(distances, loss_inh_los, 'g-', linewidth=2, label='InH LOS')
ax.plot(distances, loss_inh_nlos, 'g--', linewidth=2, label='InH NLOS')
ax.set_xlabel('Distance (m)')
ax.set_ylabel('Path Loss (dB)')
ax.set_title('3GPP Path Loss Models @ 28GHz')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# 图2: 频率对路径损耗的影响
ax = axes[0, 1]
d_fixed = 100  # m
colors = plt.cm.viridis(np.linspace(0, 1, len(frequencies)))

loss_vs_freq = []
for f in frequencies:
    loss = path_loss_3gpp_umi_los(f, d_fixed)
    loss_vs_freq.append(loss)

ax.bar(range(len(frequencies)), loss_vs_freq, color=colors, alpha=0.7)
ax.set_xticks(range(len(frequencies)))
ax.set_xticklabels(freq_labels, rotation=45, ha='right')
ax.set_ylabel('Path Loss (dB)')
ax.set_title(f'Path Loss vs Frequency @ d={d_fixed}m\n(UMi LOS)')
ax.grid(True, alpha=0.3, axis='y')

# 添加数值标签
for i, loss in enumerate(loss_vs_freq):
    ax.text(i, loss + 2, f'{loss:.0f}', ha='center', fontsize=9)

# 图3: 大气吸收和雨衰
ax = axes[1, 0]
f_range = np.linspace(1, 100, 100)  # GHz
d_km = 1  # km

# 大气吸收
atm_loss = [atmospheric_absorption(f, d_km) for f in f_range]

# 不同降雨率下的雨衰
rain_rates = [0, 5, 25, 50]  # mm/h
rain_labels = ['No Rain', 'Light (5mm/h)', 'Moderate (25mm/h)', 'Heavy (50mm/h)']
colors = ['b', 'g', 'orange', 'r']

for rain_rate, label, color in zip(rain_rates, rain_labels, colors):
    rain_loss = [rain_attenuation(f, d_km, rain_rate) for f in f_range]
    ax.plot(f_range, rain_loss, color=color, linewidth=2, label=label)

ax.plot(f_range, atm_loss, 'k--', linewidth=2, alpha=0.7, label='Atmospheric Absorption')
ax.set_xlabel('Frequency (GHz)')
ax.set_ylabel('Attenuation (dB)')
ax.set_title(f'Atmospheric and Rain Attenuation @ d={d_km}km')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_yscale('log')

# 图4: 穿透损耗比较
ax = axes[1, 1]
materials = list(penetration_loss_3gpp.keys())
lf_losses = [penetration_loss_3gpp[m]['LF'] for m in materials]
hf_losses = [penetration_loss_3gpp[m]['HF'] for m in materials]

x = np.arange(len(materials))
width = 0.35

ax.bar(x - width/2, lf_losses, width, label='Low Freq (<6GHz)', alpha=0.8)
ax.bar(x + width/2, hf_losses, width, label='High Freq (>6GHz)', alpha=0.8)

ax.set_xlabel('Material')
ax.set_ylabel('Penetration Loss (dB)')
ax.set_title('Penetration Loss: Low vs High Frequency')
ax.set_xticks(x)
ax.set_xticklabels(materials, rotation=45, ha='right')
ax.legend()
ax.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
output_path = os.path.join(output_dir, 'mmwave_propagation.png')
plt.savefig(output_path, dpi=150, bbox_inches='tight')
plt.close()
print(f"Saved: {output_path}")

print("\n关键结果:")
print(f"  28GHz, 100m处:")
print(f"    UMi LOS: {path_loss_3gpp_umi_los(28, 100):.1f} dB")
print(f"    UMi NLOS: {path_loss_3gpp_umi_nlos(28, 100):.1f} dB")
print(f"    UMa LOS: {path_loss_3gpp_uma_los(28, 100):.1f} dB")
print(f"    UMa NLOS: {path_loss_3gpp_uma_nlos(28, 100):.1f} dB")
print(f"    InH LOS: {path_loss_3gpp_inh_los(28, 100):.1f} dB")
print(f"    InH NLOS: {path_loss_3gpp_inh_nlos(28, 100):.1f} dB")

七、工程应用案例

7.1 5G网络规划案例

场景描述：某城市商业区需要部署5G基站，工作频段3.5GHz，要求覆盖半径500m，边缘速率100Mbps。

传播模型选择：

• 宏站覆盖：使用3GPP UMa模型
• 街道峡谷：使用Walfisch-Ikegami模型
• 室内深度覆盖：考虑穿透损耗

链路预算分析：

参数	数值
发射功率	46 dBm (40W)
天线增益	18 dBi
馈线损耗	3 dB
EIRP	61 dBm
路径损耗(500m, NLOS)	约135 dB
穿透损耗	15 dB
阴影衰落余量	8 dB
接收机灵敏度	-95 dBm
链路余量	约8 dB

基站布局建议：

• 站间距：300-400m
• 天线高度：25-35m
• 采用3扇区配置
• 重点区域部署微站补充覆盖

7.2 室内分布系统设计案例

场景描述：某大型购物中心(5层，每层10000m²)需要设计4G/5G室内分布系统。

设计步骤：

1. 覆盖目标确定

   • 覆盖率：>95%
   • 边缘场强：>-85dBm
   • 数据速率：下行50Mbps，上行10Mbps

2. 传播模型选择

   • 使用室内对数距离模型
   • 路径损耗指数n=2.5
   • 阴影衰落标准差σ=6dB

3. 天线布点规划

   • 吸顶天线覆盖半径：15-20m
   • 天线间距：25-30m
   • 每层天线数量：约16个

4. 功率分配设计

   • 信源功率：43dBm
   • 干放增益：30dB
   • 馈线损耗：约10dB
   • 天线口功率：约10dBm

5. 系统仿真验证

   • 使用射线追踪或主导路径模型
   • 验证覆盖均匀性
   • 检查切换区域

7.3 毫米波固定无线接入(FWA)案例

场景描述：使用28GHz毫米波为郊区住宅提供宽带接入。

系统参数：

• 基站：4T4R，EIRP 60dBm
• 用户终端：室外安装，高增益天线(24dBi)
• 目标覆盖距离：300m

传播特性分析：

• 自由空间损耗(300m, 28GHz)：109 dB
• 大气吸收：可忽略
• 雨衰余量(99.9%可用性)：5 dB
• 植被损耗：3-5 dB

链路预算：

• 总损耗：约120 dB
• 接收功率：60 + 24 - 120 = -36 dBm
• 接收机灵敏度：-80 dBm
• 链路余量：44 dB，支持高阶调制(256QAM)

部署要点：

• 视距(LOS)传播至关重要
• 需要精确的天线对准
• 考虑季节性植被变化
• 备用链路规划应对遮挡

八、总结

本教程系统介绍了城市与室内电波传播的理论基础和建模方法，包括：

1. 基础理论：自由空间传播、大尺度衰落、小尺度衰落模型
2. 城市传播：射线追踪、主导路径模型、街道峡谷效应
3. 室内传播：穿透损耗、楼层衰减、室内射线追踪
4. 毫米波特性：5G/6G高频段传播特点、波束管理
5. 工程应用：网络规划、室内分布系统、固定无线接入

关键要点：

• 不同环境需要选择合适的传播模型
• 射线追踪提供高精度但计算复杂
• 经验模型适合快速估算但精度有限
• 毫米波传播高度依赖视距条件
• 穿透损耗是室内外覆盖的关键挑战

随着5G/6G和物联网的发展，电波传播建模将继续向高精度、实时化、智能化方向发展，结合机器学习和数字孪生技术，为无线通信系统优化提供更强大的支撑。