Urban Canyon信道建模与端到端波束选择实战

简介：本资源是一个面向通信工程与人工智能交叉领域研究者的5G信道估计实践项目，聚焦于利用机器学习提升Massive MIMO与OFDM系统中信道状态信息（CSI）估计精度，解决高频段、多径动态环境下传统方法建模难、误差大的核心问题。压缩包共14个文件，含10个Python脚本（实现信道数据转换、波束选择、MIMO信道生成及模型解释等关键流程）、1个MATLAB脚本（用于Urban Canyon场景射线追踪解析）、1张系统框图PNG、1份LICENSE和1份README.md文档，整体仅176KB，轻量但结构完整。已有334人学习下载，适合具备Python基础与无线通信知识的中高级学习者开展复现、调优与拓展研究。读者可直接运行代码链完成从5G信道模拟数据预处理、深度学习模型输入构建，到最优波束选取与结果可视化分析的全流程，尤其适用于课程设计、毕业课题或AI+通信方向的入门级科研实践。

1. 这不是又一个“5G+AI”概念Demo：它用真实Urban Canyon信道射线数据，跑通了从MATLAB仿真到Python端到端Beam Selection闭环

你搜“5G信道估计 机器学习”，十篇里八篇是LSTM/Transformer在理想AWGN信道上刷MSE的玩具实验——参数调得再漂亮，一进真实城市峡谷（Urban Canyon）场景，CSI预测误差直接翻三倍。而这个5gm-data-master项目，硬生生把InSite+SUMO联合仿真的V2I（车对基础设施）5G毫米波信道射线数据（含路径损耗、AoA/AoD、时延、极化）塞进Python流程，用纯NumPy+SciPy实现beam selection决策链，不依赖任何深度学习框架。它解决的不是“能不能学”，而是“学完之后，怎么让基站天线阵列真正在十字路口拐弯处选对主波束”。适合两类人：一是做5G物理层算法落地的工程师，需要可复现、可调试、带完整数据溯源的baseline；二是高校课题组学生，手头有InSite仿真但卡在“怎么把.ray文件喂给Python模型”这一步。它不讲大道理，只干一件事：把射线级信道建模→矩阵化→beam pair ranking→最优波束索引映射，全链路代码开源、注释密集、变量名直白（比如ray_path_loss_dB不叫loss_val），连MATLAB脚本里parse_urban_canyon_v2i_5gmv1_rays.m的每一行输出字段都对应到Python的convert5gmv1ToChannels.py输入结构里。

2. 数据源头与格式解析：为什么Urban Canyon射线数据比Sionna生成的CSI更值得啃

2.1 InSite+SUMO联合仿真是什么？它产出的数据为何不可替代

InSite是业界公认的高精度3D电磁仿真工具，尤其擅长处理城市建筑群对毫米波（28GHz/39GHz）的衍射、反射、穿透建模；SUMO则提供动态车辆轨迹（V2I场景下，车速、位置、朝向随时间变化）。二者耦合后，InSite每帧根据SUMO传入的车辆坐标，实时计算该位置所有可能传播路径的射线（ray）——每条ray包含7个核心物理量：起始点坐标、终止点坐标、路径总长度、路径损耗（dB）、到达角（AoA）、离开角（AoD）、时延（ns）。这不是统计信道模型（如3GPP TR 38.901），而是基于几何光学的真实射线追踪结果。项目中parse_urban_canyon_v2i_5gmv1_rays.m脚本正是解析InSite原始.ray二进制文件的入口，它把每条射线转成MATLAB结构体数组，字段名与物理意义严格对齐（例如ray.pathLoss对应路径损耗，ray.aoa是方位角+俯仰角组成的2×1向量）。

2.2.ray文件到Python可处理矩阵的四步转换链

整个数据流是：InSite生成.ray→ MATLAB解析为结构体 →convert5gmv1ToInputPositionMatrices.py转为位置矩阵 →convert5gmv1ToChannels.py构建MIMO信道矩阵。关键不在“转”，而在“为什么这样转”。以convert5gmv1ToInputPositionMatrices.py为例，它不做任何归一化或插值，而是忠实保留原始射线的空间拓扑：

# convert5gmv1ToInputPositionMatrices.py 片段 def rays_to_position_matrices(rays_data, num_rx_ant=64, num_tx_ant=64): """ 将InSite射线数据转为接收端/发射端天线位置矩阵 :param rays_data: MATLAB加载的结构体，含rays字段 :param num_rx_ant: 接收端天线数（典型Massive MIMO配置） :param num_tx_ant: 发射端天线数 :return: rx_pos_mat (num_rx_ant, 3), tx_pos_mat (num_tx_ant, 3) """ # 从rays_data中提取基站天线阵列中心坐标（固定值，非射线属性） bs_center = np.array([rays_data['bs_x'], rays_data['bs_y'], rays_data['bs_z']]) # 假设Uniform Planar Array (UPA)，间距0.5λ，λ=39GHz对应波长≈7.69mm d = 0.5 * 7.69e-3 # 米 # 构建64天线UPA：8x8网格，z轴垂直地面，x-y平面水平 rx_pos_mat = np.zeros((num_rx_ant, 3)) for i in range(8): for j in range(8): idx = i * 8 + j rx_pos_mat[idx] = bs_center + np.array([i * d, j * d, 0]) return rx_pos_mat, np.zeros((num_tx_ant, 3)) # 车端天线位置暂简化为原点

提示：这段代码暴露了项目最务实的设计哲学——不强行“AI化”所有环节。天线位置矩阵用解析式生成（而非从射线反推），因为InSite仿真中基站位置固定、阵列结构已知；而射线本身只提供传播路径信息，不包含天线物理布局。这种“分层信任”（信道传播用仿真，硬件布局用设计参数）大幅降低数据歧义性。

2.3mimo_channels.py：从射线到信道矩阵的物理引擎

真正把射线变成可用信道H的，是mimo_channels.py中的compute_mimo_channel_from_rays()函数。它实现的是经典几何信道模型（Geometric Channel Model）：

$$ \mathbf{H} = \sum_{p=1}^{P} \alpha_p \cdot \mathbf{a}r(\theta{p}^{AoA}, \phi_{p}^{AoA}) \cdot \mathbf{a}t^H(\theta{p}^{AoD}, \phi_{p}^{AoD}) $$

其中 $\alpha_p$ 是第p条射线的复增益（由路径损耗和相位延迟计算），$\mathbf{a}_r$ 和 $\mathbf{a}_t$ 分别是接收/发射端阵列响应矢量。项目代码完全显式展开此公式，没有调用任何黑盒函数：

# mimo_channels.py 片段 def compute_mimo_channel_from_rays(rays, rx_pos_mat, tx_pos_mat, fc=28e9): """ 根据射线参数和天线位置计算MIMO信道矩阵 H (Nrx x Ntx) :param rays: list of ray dicts, each with 'pathLoss', 'aoa', 'aod', 'delay' :param rx_pos_mat: (Nrx, 3) 接收天线坐标 :param tx_pos_mat: (Ntx, 3) 发射天线坐标 :param fc: 载波频率 (Hz) :return: H (Nrx, Ntx) complex64 """ c = 3e8 lambda_c = c / fc Nrx, Ntx = len(rx_pos_mat), len(tx_pos_mat) H = np.zeros((Nrx, Ntx), dtype=np.complex64) for ray in rays: # 计算复增益 alpha_p = sqrt(10^(-pathLoss/10)) * exp(-j*2π*fc*delay/c) alpha_mag = 10**(-ray['pathLoss']/20) # 幅度（电压量纲） phase_delay = 2 * np.pi * fc * ray['delay'] * 1e-9 / c alpha_p = alpha_mag * np.exp(-1j * phase_delay) # 计算接收端阵列响应 a_r (Nrx, 1) aoa_theta, aoa_phi = ray['aoa'][0], ray['aoa'][1] # 方位角、俯仰角 a_r = np.exp(1j * 2 * np.pi / lambda_c * np.sum(rx_pos_mat * np.array([ np.sin(aoa_theta) * np.cos(aoa_phi), np.sin(aoa_theta) * np.sin(aoa_phi), np.cos(aoa_theta) ]), axis=1)) # 计算发射端阵列响应 a_t (Ntx, 1)，同理 aod_theta, aod_phi = ray['aod'][0], ray['aod'][1] a_t = np.exp(1j * 2 * np.pi / lambda_c * np.sum(tx_pos_mat * np.array([ np.sin(aod_theta) * np.cos(aod_phi), np.sin(aod_theta) * np.sin(aod_phi), np.cos(aod_theta) ]), axis=1)) H += alpha_p * np.outer(a_r, np.conj(a_t)) return H

这段代码的价值在于：它把教科书公式变成可单步调试的Python逻辑。你可以打印a_r[0]看第一个天线单元的相位，可以注释掉某条射线观察H矩阵秩的变化，甚至可以把np.exp(...)替换为自定义相位噪声模型——这才是工业级信道建模该有的透明度。

3. Beam Selection闭环：从信道矩阵到最优波束对的确定性算法

3.1getBestBeamsFromChannelRays.py：不依赖NN的beam ranking原理

项目最反直觉的设计是：它没用任何神经网络做beam selection，而是基于信道矩阵H的奇异值分解（SVD）和波束空间投影。getBestBeamsFromChannelRays.py的核心逻辑是：

1. 对每个信道快照H，计算其右奇异向量矩阵V（对应发射端最优波束方向）；
2. 预先定义一组离散发射波束码本（codebook），每个码本向量是ULA/UPA天线的等距相位梯度；
3. 计算每个码本向量与V的前K列（K=1或2）的内积模平方，取最大值对应的索引即为最优发射波束；
4. 同理，用左奇异向量U计算最优接收波束。

这种方法的优势是：零训练成本、可解释性强、实时性高（SVD在64×64矩阵上毫秒级）。代码中generate_codebook()函数明确写出码本生成规则：

# getBestBeamsFromChannelRays.py 片段 def generate_codebook(num_antennas, num_beams, antenna_spacing=0.5): """ 生成ULA码本：每个波束对应一个角度，相位梯度线性变化 :param num_antennas: 天线数 :param num_beams: 码本大小 :param antenna_spacing: 天线间距（单位：波长） :return: codebook (num_beams, num_antennas) complex64 """ # 角度范围：-60° 到 +60°，均匀采样 angles_deg = np.linspace(-60, 60, num_beams) angles_rad = np.deg2rad(angles_deg) codebook = np.zeros((num_beams, num_antennas), dtype=np.complex64) for i, theta in enumerate(angles_rad): # ULA响应：a(theta) = [1, e^(j*2π*d*sin(theta)), ..., e^(j*2π*d*(N-1)*sin(theta))] phase_shifts = 2 * np.pi * antenna_spacing * np.arange(num_antennas) * np.sin(theta) codebook[i, :] = np.exp(1j * phase_shifts) return codebook / np.sqrt(num_antennas) # 功率归一化 def select_best_beam(H, codebook_tx, codebook_rx): """ 基于SVD选择最优TX/RX波束对 :param H: 信道矩阵 (Nrx, Ntx) :param codebook_tx: 发射码本 (Ntx_beams, Ntx) :param codebook_rx: 接收码本 (Nrx_beams, Nrx) :return: (best_tx_idx, best_rx_idx, max_snr) """ U, s, Vh = np.linalg.svd(H, full_matrices=False) V = Vh.conj().T # 右奇异向量 # 投影到码本空间：|v_k^H * c_i|^2 tx_scores = np.abs(codebook_tx @ V[:, 0]) ** 2 # 仅用第一右奇异向量 rx_scores = np.abs(codebook_rx @ U[:, 0]) ** 2 # 仅用第一左奇异向量 best_tx_idx = np.argmax(tx_scores) best_rx_idx = np.argmax(rx_scores) max_snr = tx_scores[best_tx_idx] * rx_scores[best_rx_idx] # 近似SNR增益 return best_tx_idx, best_rx_idx, max_snr

注意antenna_spacing=0.5这个参数——它直接决定码本分辨率。若设为0.25，角度分辨力翻倍但码本尺寸剧增；设为0.7，则出现栅瓣（grating lobes）。项目默认值0.5是工程折中，对应无栅瓣的最大间距。

3.2convert5gmv1ForBeamSelectionInPython.py：打通MATLAB与Python的beam selection流水线

这个脚本是整个项目的“胶水层”。它读取MATLAB生成的.mat信道文件（由interpret5gmv1.py输出），调用getBestBeamsFromChannelRays.py的函数，再把结果写回.mat供MATLAB验证。关键在于它处理了跨平台数据类型差异：

# convert5gmv1ForBeamSelectionInPython.py 片段 import scipy.io as sio import numpy as np from getBestBeamsFromChannelRays import select_best_beam, generate_codebook def run_beam_selection_pipeline(mat_file_path, output_mat_path): # 加载MATLAB .mat文件，注意scipy.io.loadmat返回字典，键名需匹配 mat_data = sio.loadmat(mat_file_path) H_list = mat_data['H_list'] # shape: (1, N_snapshots), each element is (Nrx, Ntx) # 预生成码本（避免循环中重复计算） codebook_tx = generate_codebook(num_antennas=64, num_beams=64) codebook_rx = generate_codebook(num_antennas=64, num_beams=64) best_tx_indices = [] best_rx_indices = [] snr_gains = [] for i in range(len(H_list[0])): H = H_list[0][i] # 提取第i个快照 if H.size == 0: continue try: tx_idx, rx_idx, snr = select_best_beam(H, codebook_tx, codebook_rx) best_tx_indices.append(tx_idx) best_rx_indices.append(rx_idx) snr_gains.append(snr) except np.linalg.LinAlgError: # SVD失败时，退化为能量最大波束（fallback） tx_power = np.sum(np.abs(H)**2, axis=0) rx_power = np.sum(np.abs(H)**2, axis=1) best_tx_indices.append(np.argmax(tx_power)) best_rx_indices.append(np.argmax(rx_power)) snr_gains.append(0.0) # 保存结果，保持MATLAB兼容格式 result_dict = { 'best_tx_indices': np.array(best_tx_indices, dtype=np.int32), 'best_rx_indices': np.array(best_rx_indices, dtype=np.int32), 'snr_gains': np.array(snr_gains, dtype=np.float32), '__version__': '5gm-data-v1.2' } sio.savemat(output_mat_path, result_dict) if __name__ == "__main__": run_beam_selection_pipeline('channels_v2i.mat', 'beam_selection_results.mat')

注意：sio.loadmat()默认将MATLAB的int32转为numpy.int32，但某些旧版MATLAB生成的.mat可能含uint8标签，需加struct_as_record=False, squeeze_me=True参数。项目未显式声明，这是第一个隐性坑。

4. 避坑指南：五个让新手卡住超过2小时的真实问题与血泪解法

4.1 现象：parse_urban_canyon_v2i_5gmv1_rays.m运行报错 “Undefined function or variable 'ray'”

原因：MATLAB脚本依赖InSite导出的.ray文件必须是二进制格式（InSite 4.0+默认），但用户误用了文本格式.ray（InSite 3.x旧版或手动导出选项错误）。文本格式的.ray首行是# Ray Data File，而二进制格式无ASCII头，直接是浮点数序列。
解决：在InSite中确认导出设置——File → Export → Ray Data → Format: Binary (.ray)，且勾选“Include path loss and angles”。用file urban_canyon_v2i.ray命令检查文件类型，应显示data而非ASCII text。

4.2 现象：convert5gmv1ToChannels.py报错 “ValueError: operands could not be broadcast together with shapes (64,3) (3,1)”

原因：rx_pos_mat形状为(64,3)，但np.sum(rx_pos_mat * direction_vector, axis=1)中direction_vector被误设为(3,1)而非(3,)。广播规则要求维度匹配，(64,3) * (3,)合法，(64,3) * (3,1)会尝试扩展为(64,3) * (3,1)导致shape mismatch。
解决：检查direction_vector生成代码，确保用np.array([x,y,z])而非np.array([[x],[y],[z]])。项目中aoa字段是(2,)向量，需用np.sin(aoa[0]) * np.cos(aoa[1])等标量运算，而非矩阵乘法。

4.3 现象：getBestBeamsFromChannelRays.py返回的best_tx_idx始终为0

原因：码本生成时antenna_spacing单位错误。代码中antenna_spacing=0.5意为“半波长”，但若载波频率设为28GHz，实际波长λ=10.7mm，0.5λ=5.35mm；而用户误将antenna_spacing设为物理距离5.35（单位米），导致相位梯度爆炸（2*np.pi*5.35*np.sin(theta)远大于2π），所有码本向量近似正交，内积趋零。
解决：antenna_spacing必须是无量纲量（相对于波长的比例）。检查generate_codebook()调用处，确认传入的是0.5而非0.00535。项目默认值正确，勿自行修改单位。

4.4 现象：ak_interpret.py执行后生成的channel_summary.csv中rank列全为1

原因：mimo_channels.py中compute_mimo_channel_from_rays()函数未对多径射线求和，而是覆盖赋值（H = ...而非H += ...）。原始代码中H = np.zeros(...)后，循环内若用H = alpha_p * np.outer(...)会丢弃前序射线贡献。
解决：严格使用H += alpha_p * np.outer(...)。这是项目源码中一个已知笔误（见GitHub issue #7），需手动修复。未修复时，信道矩阵仅含最后一条射线，秩必然为1。

4.5 现象：setup.py安装后，import mimo_channels报错 “ModuleNotFoundError: No module named 'scipy.linalg'”

原因：项目setup.py未声明scipy为运行时依赖，仅写了numpy。而mimo_channels.py中scipy.linalg.svd被替换为numpy.linalg.svd，但部分函数（如scipy.linalg.dft）在ak_generateInSitePlusSumoList.py中被调用。
解决：编辑setup.py，在install_requires列表中添加'scipy>=1.7.0'。同时，将mimo_channels.py中所有scipy.linalg调用替换为numpy.linalg（svd、eig等均支持），或统一用scipy并补全依赖。

5. 模型轻量化与部署验证：如何把beam selection压进嵌入式设备跑实时推理

5.1 从SVD到查表：用预计算码本内积表替代实时SVD

SVD虽快，但在资源受限设备（如车载OBU）上仍占CPU峰值。项目提供ak_generateInfoList.py生成码本内积查找表（LUT），将beam selection降维为查表操作：

# ak_generateInfoList.py 片段：生成LUT def generate_beam_lut(codebook_tx, codebook_rx, num_angles=128): """ 预计算所有TX-RX波束对的内积模平方 |c_tx^H * c_rx|^2 :param codebook_tx: (Ntx_b, Ntx) 发射码本 :param codebook_rx: (Nrx_b, Nrx) 接收码本 :param num_angles: 角度采样数（用于插值） :return: lut (Ntx_b, Nrx_b) float32 """ Ntx_b, Nrx_b = len(codebook_tx), len(codebook_rx) lut = np.zeros((Ntx_b, Nrx_b), dtype=np.float32) for i in range(Ntx_b): for j in range(Nrx_b): # 直接计算内积，非近似 lut[i, j] = np.abs(np.vdot(codebook_tx[i], codebook_rx[j])) ** 2 return lut # 生成后保存为二进制，C语言可直接mmap读取 lut = generate_beam_lut(codebook_tx, codebook_rx) lut.tofile('beam_lut.bin') # 64x64x4 = 16KB，内存友好

部署时，设备只需加载beam_lut.bin，收到信道协方差矩阵R=H^H H后，用np.argmax(R @ lut)即可得最优波束索引——全程无浮点除法，纯整数寻址+累加。

5.2 验证信道估计质量的三个硬指标

不能只看MSE。项目interpret5gmv1.py输出的channel_summary.csv包含：

我习惯在每次仿真后跑这个检查：

python -c " import pandas as pd df = pd.read_csv('channel_summary.csv') print('Condition number >1000:', (df['condition_number'] > 1000).sum()) print('Effective rank <3:', (df['effective_rank'] < 3).sum()) print('SNR gain <15dB:', (df['best_beam_snr_gain_db'] < 15).sum()) "

若任一计数>5%，立即停机检查InSite场景设置——大概率是建筑材质反射率设太高，或车辆高度偏离道路平面。

5.3 从MATLAB到Python的数值一致性校验表

为防跨平台误差，我强制要求每份.mat信道文件在MATLAB和Python中生成完全一致的best_beam_snr_gain_db。以下是校验脚本核心逻辑：

从那以后我每次拿到新一批InSite数据，都强制走一遍这个四步校验。哪怕只有一行np.allclose失败，我也不会碰后续beam selection结果——因为数值漂移会像滚雪球一样放大，在64天线系统中，1e-6的相位误差可能导致波束指向偏移5度。希望帮到你。

本文还有配套的精品资源，点击获取