波束成形、DOA估计与RIS联合仿真:从阵列模型到算法实现
2026/9/15 5:24:57 网站建设 项目流程

简介:这是一份基于MATLAB的波束成形、DOA估计与RIS(可重构智能表面)技术仿真资源,面向通信、电子信息、数学等专业学生的课程设计、期末大作业与毕业设计场景。代码采用参数化编程,参数可灵活调整,注释清晰,并附带可直接运行的案例数据,能帮助使用者快速掌握阵列信号处理与智能超表面辅助通信的核心原理。压缩包共19个文件,以18个.m源码文件为主,另含1个.md说明文档,整体仅39KB,轻量精简、便于本地部署与二次开发。目前已有92人学习或下载,适合正在学习相关课程或计划开展科研复现的学生参考:通过阅读代码与案例,可系统理解波束成形、DOA估计及RIS辅助通信等关键步骤,也可在此基础上修改参数、替换场景,完成自己的仿真实验。

1. 一个 ZIP 里同时出现三个词,先别急着解压

拿到“波束成形-DOA-RIS.zip”这个包,多数人的第一反应是解压、找 .m 或 .py 文件直接运行。但真正决定这个包价值的是三个英文缩写之间的依存关系:波束成形决定能量怎么打出去,DOA 估计决定朝哪个方向打,RIS 则是中途用来改变电磁传播路径的一张“可编程反射镜”。三者放在同一个压缩包里,通常不是三个独立 demo 的简单合并,而是一条“先探测信道 → 估计用户方位 → 设计收发权值 → 协同调 RIS 相位”的链路。这篇文章按我拿到此类仿真包时会走的路径来讲:先固定阵列模型,再做 DOA,再引入 RIS 联合设计,最后落到验证方法。无论包里是 MATLAB 还是 Python 工程,读完之后你至少能判断这段代码在仿什么、参数改哪里、结果是否可信。

2. 阵列模型与波束成形:把坐标系和导向矢量先钉死

2.1 均匀线阵的导向矢量:窄带、远场、平面波是默认前提

在波束成形和 DOA 估计的代码里出现最多的结构是均匀线阵(ULA)。设阵元数为 (M),阵元间距为 (d),信号入射角为 (\theta),采用窄带、远场、平面波假设时,第 (m) 个阵元相对参考阵元的传播时延为 (m d \sin\theta / c)。用载波频率归一化后,导向矢量写成:

[ \mathbf{a}(\theta) = \left[1,\ e^{j 2\pi (d/\lambda) \sin\theta},\ \ldots,\ e^{j 2\pi (M-1)(d/\lambda) \sin\theta}\right]^T ]

(d/\lambda) 是最关键的归一化参数。取 0.5(半波长)时阵列无栅瓣,同时角度分辨率达到常用配置的平衡点。(d/\lambda) 小于 0.4 时主瓣变宽;大于 0.5 时视场边缘出现栅瓣,DOA 估计会出现模糊峰。波束成形与 DOA 这两个模块的精度上限,最早就由这个比值决定。

坐标系约定同样值得花三十秒确认。多数代码以阵列法线为 0°,一侧为正、另一侧为负;少数实现用端射方向做参考。若把 (\theta) 直接代入 (\sin\theta),两侧符号换掉,波束图会左右翻转,MUSIC 谱的峰也会出现在错误象限。拿到包的第一步,先找 angle_grid 的定义,确认角度范围和步进——这个动作能避免后面几十分钟的无效调试。

2.2 从匹配滤波器到 MVDR:两种权值设计的取舍

常规相控阵直接取 (\mathbf{w} = \mathbf{a}(\theta_0)),相当于对目标方向的信号做相干累加,实现简单但对干扰没有任何抑制能力。MVDR(最小方差无失真响应)在目标方向保持增益 1,同时最小化输出总功率,从而把干扰方向“压出凹口”。其闭式解为:

[ \mathbf{w} = \frac{\mathbf{R}^{-1} \mathbf{a}(\theta_0)}{\mathbf{a}(\theta_0)^H \mathbf{R}^{-1} \mathbf{a}(\theta_0)} ]

这里 (\mathbf{R}) 是阵列接收数据的协方差矩阵。实际工程中快拍数不足时样本协方差矩阵求逆会放大噪声,通常要加对角加载:(\mathbf{R} + \lambda_0 \mathbf{I})。加载系数 (\lambda_0) 一般取 (\mathbf{R}) 对角线均值乘以 (10^{-2}) 到 (10^{-1})。过小失去稳定作用,过大会退化成匹配滤波器。代码包里如果只有 (\mathbf{a}(\theta_0)) 直接加权而没有矩阵求逆,那就是常规相控阵,谈不上自适应抗干扰。

2.3 用 Python 跑通最小波束成形闭环并检查主瓣指向

import numpy as np def steering_vector(M, theta_deg, d_over_lambda=0.5): """均匀线阵导向矢量,theta_deg 是相对阵列法线的角度""" theta = np.deg2rad(theta_deg) m = np.arange(M) return np.exp(1j * 2 * np.pi * d_over_lambda * m * np.sin(theta)) def mvdr_weights(Rxx, a): """MVDR 权值,带对角加载防止 R 奇异""" lam = 1e-2 * np.trace(Rxx) / Rxx.shape[0] Rinv = np.linalg.inv(Rxx + lam * np.eye(Rxx.shape[0])) return Rinv @ a / (a.conj().T @ Rinv @ a) M, N = 8, 500 theta0, snr_db = 10, 10 s = (np.random.randn(N) + 1j*np.random.randn(N)) / np.sqrt(2) a0 = steering_vector(M, theta0) noise_power = 10**(-snr_db/10) X = np.outer(a0, s) + np.sqrt(noise_power) * ( np.random.randn(M, N) + 1j*np.random.randn(M, N)) / np.sqrt(2) Rxx = (X @ X.conj().T) / N w = mvdr_weights(Rxx, a0) angles = np.linspace(-90, 90, 361) pattern = [abs(w.conj() @ steering_vector(M, ang))**2 for ang in angles] peak = angles[np.argmax(pattern)] print(f"主瓣峰值: {peak}°, 期望: {theta0}°")

逻辑说明:先按“单信号加高斯白噪声”生成 8 阵元、500 快拍的数据,样本协方差 Rxx 是 (8\times8) 的复数矩阵。MVDR 部分做了一次对角加载,加载系数取对角线均值的 1%,对 10 dB 信噪比场景足够稳定。归一化由复数除法完成,确保期望方向增益为 1。

参数可以按场景来改:阵元数 M 决定空间自由度,M 越大旁瓣越低;N 是快拍数,从 500 降到 50 时,Rxx 估计误差变大,主瓣指向偏移会更明显,这也是信道估计类项目里最常见的坑。把 snr_db 改成 -5,再看主瓣峰值位置,能直观理解“低信噪比时波束成形精度不再由孔径决定,而由协方差估计误差决定”。

3. DOA 估计与 MUSIC 的实用参数:子空间维度是第一个检查点

3.1 为什么 MUSIC 能做出超分辨率

MUSIC 类方法的基本前提是信号子空间与噪声子空间正交。对阵列协方差 (\mathbf{R}) 做特征分解后,大特征值对应的特征向量张成信号子空间,余下特征向量张成噪声子空间 (\mathbf{E}_n)。由于导向矢量落在信号子空间内,它与噪声子空间正交,因此在真实入射角附近满足 (\mathbf{a}(\theta)^H \mathbf{E}_n \approx \mathbf{0})。空间谱定义为该内积模长的倒数:

[ P_{\text{MUSIC}}(\theta) = \frac{1}{\mathbf{a}(\theta)^H \mathbf{E}_n \mathbf{E}_n^H \mathbf{a}(\theta)} ]

谱峰位置即 DOA 估计结果。超分辨能力的来源在于:特征分解能区分角度差小于阵列瑞利波束宽度的两个信号源,只要快拍数和信噪比足够。MUSIC 谱不是功率谱,它没有物理单位,谱峰高度没有绝对意义,只有相对意义。

注意一个微妙点:MUSIC 能分辨多少个源,取决于样本协方差矩阵能支持多少个维度,而不是真实存在多少个来波。信源数 K 估多了,本应属于噪声子空间的向量被划入信号子空间,谱会漏峰;K 估少了,伪峰出现。常见做法是用 MDL 或 AIC 准则估计信源数,再传给 MUSIC。代码包里如果直接写死 K=1,而仿真场景里有两个用户,结果必然失真。

3.2 SubspaceNet 的思路:用网络拟合子空间投影

如果做 5G-A 或通感一体化这类场景,快拍数经常少于阵元数,此时样本协方差不满秩,特征分解的性能边界迅速退化。SubspaceNet 这类方法的思路是用深度网络直接学习从阵列数据到信号子空间投影的映射,输入可以是复协方差矩阵的实虚部堆叠,监督信号则由高快拍数下的特征分解产生。训练完成后推理时避开显式特征分解,在低快拍时比传统 MUSIC 更稳。放在 ZIP 包里的位置通常是 DOA 估计模块的一个可选替换:传统代码里是 np.linalg.eigh,替换后是一个网络推理函数。实际使用中要小心训练场景与测试场景的天线结构不一致——子空间形状和阵元布局强绑定,换阵列就要重新训练或微调。

3.3 一个不依赖外部库的 MUSIC 最小实现

def music_spectrum(X, K, angles): """X: M×N 阵列数据; K: 信源个数; angles: 扫描角度列表""" M = X.shape[0] Rxx = (X @ X.conj().T) / X.shape[1] eigvals, eigvecs = np.linalg.eigh(Rxx) En = eigvecs[:, :M-K] # 噪声子空间:特征值升序排列,取前 M-K 列 spec = [] for ang in angles: a = steering_vector(M, ang) spec.append(1.0 / np.abs(a.conj() @ En @ En.conj().T @ a)) return np.array(spec) angles = np.linspace(-90, 90, 361) spec = music_spectrum(X, K=1, angles=angles) est = angles[np.argmax(spec)] print(f"DOA 估计结果: {est}°, 真实来波: {theta0}°")

逻辑说明:np.linalg.eigh 返回的特征值按升序排列,因此噪声子空间取前 M-K 列。空间谱每个角度独立计算一次倒数,扫描 361 个角度时开销不大;工程实现里可以用矩阵一次性完成计算。K 的选取直接影响谱峰质量,可以在 -10 dB 信噪比下把 K 改成 2,观察伪峰是否出现。把 M、快拍 N、信噪比三个参数串起来看会更直观:N 小于 M 时样本协方差秩亏,特征分解会把噪声能量摊到后 M-K 个特征值上,谱峰抖动明显。工程上我会先保证 N ≥ 4M,再谈超分辨率。

4. RIS 加入后模型多了一层反射:级联信道与联合设计

4.1 级联信道模型和反射系数的影响

RIS 由 N 个反射单元组成,每个单元相当于一个可调相移器,反射系数写为 (\beta_n e^{j\varphi_n}),(\beta_n) 为幅度(理想情况取 1),(\varphi_n) 为相移。在远场窄带假设下,基站到用户的等效信道由直射径 (\mathbf{h}d)、基站-RIS 信道 (\mathbf{G})((N\times M) 矩阵)、RIS-用户信道 (\mathbf{h}{ru})(N 维向量)三部分组成,级联信道表达式为:

[ \mathbf{h}_{\text{eff}} = \mathbf{h}d + \mathbf{h}{ru}^H \boldsymbol{\Phi} \mathbf{G} ]

其中 (\boldsymbol{\Phi} = \mathrm{diag}(\beta_1 e^{j\varphi_1}, \ldots, \beta_N e^{j\varphi_N}))。RIS 本身不放大信号,它通过调整相位使来自不同反射单元的径在接收端相干叠加,形成类似波束成形的“空间聚焦”。这也是为什么 RIS 相关工作里频繁出现“被动波束成形”这个词——它不消耗发射功率,却改变信道增益。

DOA 与 RIS 的结合点在于:要对用户做定向波束成形,先要知道用户方向;当直射径被遮挡时,波达方向信息只能通过 RIS 反射链路获得。因此这类代码包里会先固定一组相移做探测,用 DOA 估计出反射来波方向,再据此优化下行权值。这里的 DOA 估计对象是“经过 RIS 反射后的虚拟路径方向”,不是直接来波方向。

4.2 相移设计和两种常见联合优化策略

理想相移是无约束最优解:(\varphi_n = \arg(\mathbf{h}_{ru,n}) - \arg([\mathbf{G}\mathbf{a}(\theta)]_n)),即把每个反射单元的用户信道相位与到达相位对齐。但每个单元相移是量化的,例如 2-bit 移相时 (\varphi_n \in {0, \pi/2, \pi, 3\pi/2});量化误差会让波束成形增益下降,1-bit 量化损失约 3.9 dB,2-bit 损失约 0.9 dB。代码包中看到相位取整操作,就是这种离散化。下面给出理想相移计算并对齐的片段:

N_ris, M = 32, 8 theta_user = 30 G = (np.random.randn(N_ris, M) + 1j*np.random.randn(N_ris, M)) / np.sqrt(2) h_ru = (np.random.randn(N_ris) + 1j*np.random.randn(N_ris)) / np.sqrt(2) g_beamformed = G @ steering_vector(M, theta_user) target_phase = np.angle(h_ru * g_beamformed) phi_bits = 2 phi_levels = 2**phi_bits step = 2*np.pi / phi_levels phi_quant = np.round(target_phase / step) * step Phi = np.diag(np.exp(1j * phi_quant)) cascaded = np.abs(h_ru.conj() @ Phi @ G @ steering_vector(M, theta_user)) perfect = np.abs(h_ru.conj() @ np.diag(np.exp(1j * target_phase)) @ G @ steering_vector(M, theta_user)) print(f"量化后级联增益: {cascaded:.3f}, 理想增益: {perfect:.3f}")

逻辑说明:先计算每个反射单元的到达复增益 g_beamformed,再与用户信道 h_ru 逐元素相乘取相位,得到理想对齐相位。量化到 2-bit 后构造对角矩阵,计算级联信道增益。运行结果中量化增益会略低于理想增益,这个差值就是离散化损失。实际系统中更大的损失来自信道估计误差:当 G 和 h_ru 存在 10% 幅度误差和 5° 相位误差时,算出的相移已偏离最优值,此时继续提高相移分辨率对增益提升有限。

联合设计方面,文献里最常用的两种做法是两阶段法和交替优化。两阶段法先关闭 RIS 或使用正交导频做信道探测,完成 DOA/信道估计后一次性设计权值和相移,实现简单但存在误差传播;交替优化把问题拆成“固定 Φ 优化 w”和“固定 w 优化 Φ”两个子问题循环求解,性能更好但收敛性需要验证。zip 包里若存在 while 循环或 for iter 结构,多半是交替优化;两阶段法通常只有顺序执行的主脚本,没有迭代逻辑。

4.3 常见仿真参数设置和选择理由

以下表格是 RIS 辅助通信仿真中常见的参数范围,便于对照 zip 包里的默认配置判断是否合理:

参数常见取值说明
RIS 单元数 N16 ~ 128越大增益越高,但 CSI 开销和优化复杂度随之上涨
单元间距(\lambda/2)与天线类似,过大会出现反射栅瓣
相移量化位 b1 ~ 3 bit1-bit 增益损失明显,2-bit 是性能和复杂度平衡点
基站天线数 M8 ~ 64决定 DOA/波束成形的空间分辨率
链路信噪比-5 ~ 15 dB低于 -5 dB 时 DOA 误差对波束增益的影响显著
快拍数 N_snap128 ~ 1024匹配 DOA 高分辨需求,N_snap ≥ 4M
用户移动性低速/准静态RIS 相移更新周期一般按信道相干时间设计

表格看完后要核对 zip 包中的三个尺度关系:载频决定波长,进而决定阵元间距;信道相干时间决定相移更新间隔;DOA 估计间隔决定波束更新速率。三者若同处一个数量级,系统存在严重的闭环风险——估计完 DOA,用户已经移动,波束指向已经过时。

5. 解开 ZIP 包后,用三组指标验证链路真正闭环

5.1 先做遮挡实验:直射径为零时 DOA 是否依然可用

遮挡场景是 RIS 参与链路的最大价值所在。手动把直射增益置零,只保留反射径,跑一遍 DOA,观察谱峰是否还出现在真实入射角附近。如果谱峰消失或偏移超过 5°,说明代码包中的 DOA 模块默认依赖直射径,存在隐藏的链路假设。这一步比任何调试都值得先做,它能暴露多数代码包“跑得通但没闭环”的问题。

5.2 三组指标:估计误差、波束增益、迭代收敛

第一组指标是 DOA 估计误差,做 200 次蒙特卡洛仿真,统计估计角度的均值和标准差。在快拍数 256、信噪比 10 dB 的配置下,8 阵元 ULA 对单目标的理论 CRB 大约在 0.3° ~ 1° 量级,实际估计的标准差不应超过该数值的 3 ~ 5 倍。第二组指标是波束成形增益与理想值的差距,做法是把 MVDR 权值替换成纯匹配滤波权值,比较输出信噪比差异;若 MVDR 与匹配滤波增益差小于 1 dB,说明协方差估计不准,需要增加快拍或对角加载系数。第三组指标是 RIS 迭代收敛曲线,输出目标函数值,检查算法是否能在 50 次迭代内稳定收敛;若振荡幅度不下降,排查步长设置和相位初始化方式。

5.3 用一张图同时验证三个模块

把真实入射角、MUSIC 谱峰位置、MVDR 方向图主瓣指向画在同一张图上。正常闭环时这三个角度的偏差应在 1° 以内;出现系统性偏移(三个角度输出一致但整体偏移)说明角度定义有偏差;出现 MUSIC 峰在 0° 附近、MVDR 主瓣在另一侧的情况,说明两个模块用了不同的角度坐标系假设。这时对比代码中 angle_grid 的生成方式和 steering_vector 的符号约定,通常能快速定位问题。对代码包而言,靠人工阅读难以发现的坐标系错位,用一张图往往能立即看到——这也是我在解开任何波束成形-DOA-RIS 代码包后必做的第一件事。

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询