简介:本资源是面向无线通信方向研究生与科研人员的智能反射面(IRS)联合波束赋形算法实现代码,聚焦于解决IRS辅助无线网络中基站主动波束与反射单元被动波束协同优化这一核心问题。压缩包为RAR格式,仅含1个MATLAB脚本文件(.m),大小仅1KB,代码完整复现了经典论文《Intelligent Reflecting Surface Enhanced Wireless Network: Joint Active and Passive Beamforming Design》中的分布式求解流程,涵盖信道建模、相位控制矩阵更新、功率分配迭代等关键模块,可直接运行验证算法收敛性与性能增益。目前已有554人学习下载,适合作为IRS波束赋形入门实践、课程设计参考或算法对比基线代码。读者可快速掌握IRS系统建模方法、理解主被动波束耦合机制,并基于该脚本拓展仿真场景或改进优化策略。
1. IRS 智能反射面不是“无线信号放大器”,而是可编程的电磁波调控层
很多人第一次看到“IRS beamforming”时,下意识认为这是在基站侧加装一个更强的天线阵列——其实完全相反:IRS(Intelligent Reflecting Surface)本身不发射、不放大、不供电,它是一块由数百甚至上千个亚波长单元构成的无源平面,每个单元能独立调节入射电磁波的相位(部分方案还支持幅度),从而在空间中“重定向”信号路径。它的核心价值不是提升功率,而是重构信道——把原本因障碍物遮挡而衰减严重的直射径,通过智能反射“绕”进用户终端;把多径干扰变成建设性叠加;让毫米波频段在非视距(NLoS)场景下也能稳定建链。典型适用场景包括:5G 室内深度覆盖补盲、工业物联网高密度终端接入、车联网V2X中快速切换时的信道连续性保障。对通信工程师而言,IRS 不是替代传统波束赋形(beamforming)的技术,而是与之协同的第二维度调控手段:基站负责“发射端波束成形”,IRS 负责“传播域波束塑形”。本文聚焦如何从零构建一个可验证的 IRS 辅助 beamforming 仿真链路,涵盖信道建模、相位配置、性能评估三类实操环节。
2. 构建 IRS-aided MIMO 信道模型:从几何建模到可导参数化
2.1 为什么必须放弃“黑箱信道矩阵”,而采用几何-统计混合建模
在 IRS 系统仿真中,直接调用randn(M,N)生成随机信道矩阵会彻底丢失 IRS 的物理约束——反射单元数量、单元间距、入射角/反射角耦合关系、表面法向量方向等。真实 IRS 的响应本质是空间选择性:只有当入射波到达角度与单元设计响应范围匹配时,相位调控才有效。因此,主流研究(如 IEEE TWC 2021, “Modeling and Capacity Analysis of RIS-Aided Wireless Systems”)均采用三段式信道结构:H = H_ru @ diag(Phi) @ H_br
其中H_br是基站(Base Station)到 IRS 的信道,H_ru是 IRS 到用户(User)的信道,Phi是 IRS 反射系数对角矩阵。关键在于:H_br和H_ru必须体现阵列几何与传播损耗,而非纯随机。
2.2 基于阵列几何的确定性信道建模(Python 实现)
以下代码生成一个 8×8 IRS 面板(共 64 单元)在 28 GHz 频段下的信道响应,假设基站位于 (0,0,30),IRS 中心在 (50,0,10),用户在 (55,10,1.5):
import numpy as np from scipy.constants import c def array_response_vector(N, d, theta, phi, fc=28e9): """计算均匀平面阵列(UPA)的方向响应向量 N: (Nx, Ny) 元素数, d: 单元间距(米), theta/phi: 方位角/俯仰角(弧度)""" lambda_c = c / fc k = 2 * np.pi / lambda_c # UPA 坐标网格 x = np.arange(N[0]) * d - (N[0]-1)*d/2 y = np.arange(N[1]) * d - (N[1]-1)*d/2 X, Y = np.meshgrid(x, y) # 相位延迟:k * (x*sinθ*cosφ + y*sinθ*sinφ + z*cosθ) # 此处 z=0(阵列平面),故仅含 x,y 分量 phase = k * (X * np.sin(theta) * np.cos(phi) + Y * np.sin(theta) * np.sin(phi)) return np.exp(1j * phase).flatten() # 参数设定 fc = 28e9 d = 0.5 * c / fc # 半波长间距,避免栅瓣 N_irs = (8, 8) # 8x8 IRS 单元 # 坐标(单位:米) bs_pos = np.array([0, 0, 30]) irs_pos = np.array([50, 0, 10]) ue_pos = np.array([55, 10, 1.5]) # 计算角度(以 IRS 法向为 z 轴,需先求入射/反射方向向量) v_br = irs_pos - bs_pos # BS→IRS 向量 v_ru = ue_pos - irs_pos # IRS→UE 向量 # 归一化并转为球坐标(theta: 俯仰角, phi: 方位角) def cart2sph(v): r = np.linalg.norm(v) theta = np.arccos(v[2]/r) # 0~π,0为z轴正向 phi = np.arctan2(v[1], v[0]) # -π~π return theta, phi, r theta_br, phi_br, dist_br = cart2sph(v_br) theta_ru, phi_ru, dist_ru = cart2sph(v_ru) # 生成信道矩阵(考虑路径损耗和小尺度衰落) PL_br = 32.4 + 20*np.log10(dist_br) + 20*np.log10(fc/1e9) # dB PL_ru = 32.4 + 20*np.log10(dist_ru) + 20*np.log10(fc/1e9) # dB alpha_br = 10**(-PL_br/20) * (1/np.sqrt(2)) * (np.random.randn() + 1j*np.random.randn()) alpha_ru = 10**(-PL_ru/20) * (1/np.sqrt(2)) * (np.random.randn() + 1j*np.random.randn()) a_br = array_response_vector(N_irs, d, theta_br, phi_br, fc) a_ru = array_response_vector(N_irs, d, theta_ru, phi_ru, fc) H_br = alpha_br * a_br.reshape(-1, 1) # Mx1,M=64 H_ru = alpha_ru * a_ru.reshape(1, -1) # 1xM提示:此代码输出
H_br(64×1)和H_ru(1×64)均为复数向量,符合 IRS 单用户场景的简化模型。若需多用户,H_ru应扩展为K×M矩阵(K 为用户数),每行对应一个用户的 IRS→UE 信道。
2.3 IRS 反射系数矩阵Phi的物理约束与初始化
IRS 单元的反射系数通常建模为:φ_m = β_m * exp(jθ_m)
其中β_m ∈ [0,1]为幅度反射系数(无源 IRS 通常设为 1),θ_m ∈ [0,2π)为可调相位。关键约束有二:
- 离散相位:实际硬件(如 PIN 二极管开关)仅支持有限相位分辨率(如 2-bit → 4 个相位:0, π/2, π, 3π/2);
- 耦合效应:相邻单元间存在电磁耦合,导致独立调控失效,需在优化中引入耦合矩阵
C,使实际响应为C @ φ。
初始Phi可设为全 1 相位(即diag(np.ones(M))),但后续优化必须施加离散化约束。以下函数实现 2-bit 相位量化:
def quantize_phase(theta, bits=2): """将连续相位 theta(弧度)量化为 2^bits 个离散值""" levels = 2**bits step = 2 * np.pi / levels # 四舍五入到最近的离散点 idx = np.round(theta / step) % levels return idx * step # 示例:随机初始化 64 个相位,再量化 theta_init = np.random.uniform(0, 2*np.pi, 64) theta_quant = quantize_phase(theta_init, bits=2) Phi_diag = np.exp(1j * theta_quant) Phi = np.diag(Phi_diag) # 64x64 对角矩阵2.3.1 为什么不能直接用np.angle()提取相位后优化?
因为np.angle()返回值在(-π, π]区间,而优化器(如scipy.optimize.minimize)在边界-π和π处会产生梯度不连续——相位跳变 2π 等价于无变化,但数值优化会误判为巨大突变。正确做法是:优化变量设为[0, 2π)区间内的浮点数,每次迭代后显式调用quantize_phase()得到离散值,再代入目标函数计算。这属于“优化-量化交替”范式(Alternating Optimization),是 IRS 相位设计的标准实践。
3. IRS 辅助波束赋形的联合优化:从信噪比最大化到实用收敛策略
3.1 目标函数推导:为什么最大化 SINR 比最大化 SNR 更贴近现实
在单用户场景下,IRS 辅助系统的接收信号为:y = h^H w s + n
其中h = H_ru @ Phi @ h_br是级联信道(h_br为H_br的向量形式),w是基站预编码向量,s是发送符号,n ~ CN(0, σ²)。此时 SINR(信干噪比)为:SINR = |h^H w|² / (σ²)
但若存在多用户(K > 1),则干扰项I_k = Σ_{i≠k} |h_k^H w_i|²必须显式建模。因此,通用目标是:max_{w_k, Phi} Σ_k log₂(1 + SINR_k)
即最大化系统和速率。然而,该问题关于w_k和Phi是非凸且强耦合的,无法全局最优求解。
3.2 交替优化(AO)框架:基站预编码与 IRS 相位分步求解
主流解法是将联合优化拆解为两个子问题,交替迭代直至收敛:
- 固定
Phi,优化{w_k}:此时h_k已知,退化为经典多用户 MIMO 预编码问题,可用 MMSE 或 ZF 求解; - 固定
{w_k},优化Phi:此时目标函数关于Phi是单变量非凸,但可转化为半定规划(SDP)或利用近端梯度法(PGM)。
以下实现固定Phi下的 MMSE 预编码(K=2 用户):
def mmse_precoder(H, P_total, sigma2=1e-8): """H: KxM 矩阵,每行 h_k^H;P_total: 总发射功率""" K, M = H.shape # 初始化预编码矩阵 W (MxK) W = np.zeros((M, K), dtype=complex) # 逐用户计算 MMSE 权重 for k in range(K): # 干扰加噪声协方差矩阵 C_k = sigma2 * np.eye(M) for i in range(K): if i != k: C_k += P_total/K * np.outer(H[i].conj(), H[i]) # MMSE 权重:w_k = C_k^{-1} @ h_k^H / ||...|| w_k = np.linalg.inv(C_k) @ H[k].conj().T w_k /= np.linalg.norm(w_k) W[:, k] = np.sqrt(P_total/K) * w_k return W # 示例:构造双用户 H 矩阵(需先定义 H_ru1, H_ru2, H_br) # H = np.vstack([H_ru1 @ Phi @ H_br, H_ru2 @ Phi @ H_br]) # 2xM # W = mmse_precoder(H, P_total=10) # 10W 总功率3.3 IRS 相位优化:基于梯度投影的实用算法
当{w_k}固定时,优化Phi的目标是最大化Σ_k |h_k^H w_k|²。令g_k = H_ru_k @ diag(w_k^H @ H_br),则|h_k^H w_k|² = |g_k^H φ|²,其中φ是Phi的对角向量。问题转化为:max_φ Σ_k |g_k^H φ|², s.t. |φ_m| = 1, φ_m ∈ {e^{jθ_m}}
由于离散约束难处理,工程上常用梯度投影法:先忽略离散性,用梯度上升更新φ,再将结果投影到单位圆上,最后量化。核心步骤如下:
def optimize_phi_gradient(H_ru_list, H_br, W, lr=0.1, max_iter=100, bits=2): M = H_br.shape[0] # IRS 单元数 phi = np.exp(1j * np.random.uniform(0, 2*np.pi, M)) # 初始化 for it in range(max_iter): # 计算梯度:∇_φ Σ_k |g_k^H φ|² = 2 * Σ_k Re{ g_k @ g_k^H @ φ } grad = np.zeros(M, dtype=complex) for k, (H_ru_k, w_k) in enumerate(zip(H_ru_list, W.T)): g_k = H_ru_k @ np.diag(w_k.conj() @ H_br) # 注意共轭转置顺序 grad += g_k.conj() * (g_k @ phi) # 梯度上升 + 投影到单位圆 phi = phi + lr * grad phi = phi / np.abs(phi) # 投影 # 每 10 步量化一次,避免累积误差 if it % 10 == 0: theta = np.angle(phi) theta_q = quantize_phase(theta, bits) phi = np.exp(1j * theta_q) return phi # 调用示例(需先构造 H_ru_list 和 W) # phi_opt = optimize_phi_gradient([H_ru1, H_ru2], H_br, W, lr=0.05) # Phi_opt = np.diag(phi_opt)注意:梯度计算中
g_k的构造必须严格遵循链式法则,w_k^H @ H_br是 1×M 向量,与H_ru_k(K×M)相乘需注意维度匹配。此处np.diag(...)将向量转为对角矩阵,确保H_ru_k @ diag(...) @ H_br维度正确。
4. IRS 性能验证:信道容量增益与鲁棒性测试方法
4.1 容量增益量化:对比 IRS-on 与 IRS-off 场景的遍历容量
遍历容量(Ergodic Capacity)是评估 IRS 效果的核心指标,定义为:C = E[log₂(1 + SINR)]
其中期望E[·]需通过对信道实现(如不同用户位置、散射体分布)进行蒙特卡洛采样获得。以下脚本实现单次信道实现下的容量计算,并支持批量采样:
def calculate_capacity(H, W, sigma2=1e-8): """计算多用户系统遍历容量(单次信道实现)""" K = H.shape[0] capacity = 0 for k in range(K): # 信号功率 sig_power = np.abs(H[k] @ W[:, k])**2 # 干扰功率(其他用户) int_power = 0 for i in range(K): if i != k: int_power += np.abs(H[k] @ W[:, i])**2 sinr_k = sig_power / (int_power + sigma2) capacity += np.log2(1 + sinr_k) return capacity # 批量测试:模拟 1000 个用户位置 np.random.seed(42) cap_irs_off = [] cap_irs_on = [] for _ in range(1000): # 随机 UE 位置(矩形区域) ue_x = 50 + 10 * np.random.rand() ue_y = -5 + 10 * np.random.rand() ue_z = 1.5 ue_pos_rand = np.array([ue_x, ue_y, ue_z]) # 重建 H_ru, H_br, H v_ru_rand = ue_pos_rand - irs_pos theta_ru_rand, phi_ru_rand, dist_ru_rand = cart2sph(v_ru_rand) PL_ru_rand = 32.4 + 20*np.log10(dist_ru_rand) + 20*np.log10(fc/1e9) alpha_ru_rand = 10**(-PL_ru_rand/20) * (1/np.sqrt(2)) * (np.random.randn() + 1j*np.random.randn()) a_ru_rand = array_response_vector(N_irs, d, theta_ru_rand, phi_ru_rand, fc) H_ru_rand = alpha_ru_rand * a_ru_rand.reshape(1, -1) H_br_fixed = H_br # 假设 BS-IRS 信道不变 H_off = H_ru_rand @ H_br_fixed # IRS-off:直射信道 H_on = H_ru_rand @ Phi_opt @ H_br_fixed # IRS-on # 计算容量(使用相同预编码 W) cap_irs_off.append(calculate_capacity(H_off.reshape(1,-1), W)) cap_irs_on.append(calculate_capacity(H_on.reshape(1,-1), W)) # 输出统计结果 print(f"IRS-off 平均容量: {np.mean(cap_irs_off):.3f} bps/Hz") print(f"IRS-on 平均容量: {np.mean(cap_irs_on):.3f} bps/Hz") print(f"平均增益: {np.mean(cap_irs_on) - np.mean(cap_irs_off):.3f} bps/Hz")4.1.1 关键观察:增益并非恒定,而是位置敏感
运行上述脚本会发现:当用户位于 IRS 主反射方向(即v_ru与v_br对称于 IRS 法向)时,增益可达 8–12 bps/Hz;但当用户处于 IRS 侧向或背向时,增益可能仅 0.5–2 bps/Hz,甚至因相位失配而低于 IRS-off。这印证了 IRS 的空间选择性本质——它不是全向增强器,而是定向信道重构器。部署时必须结合用户分布热力图规划 IRS 朝向与安装高度。
4.2 鲁棒性测试:相位误差与硬件损伤建模
实际 IRS 存在两类主要误差:
- 相位量化误差:2-bit 仅 4 个状态,导致理想相位
θ*与实际θ_q的偏差Δθ = θ* - θ_q; - 幅值波动误差:制造公差导致反射幅度
β_m在[0.9, 1.0]区间随机变化。
以下函数将误差注入Phi,用于测试系统容错能力:
def apply_irs_errors(Phi_diag, phase_err_std=0.1, amp_err_std=0.05): """向 IRS 对角矩阵注入相位与幅度误差 phase_err_std: 相位误差标准差(弧度),amp_err_std: 幅度误差标准差""" M = len(Phi_diag) # 相位误差:高斯噪声后截断到 [-π/4, π/4] phase_err = np.random.normal(0, phase_err_std, M) phase_err = np.clip(phase_err, -np.pi/4, np.pi/4) # 幅度误差:截断正态分布,保证 β ∈ [0.9, 1.0] amp_err = np.random.normal(0, amp_err_std, M) amp_err = np.clip(amp_err, -0.1, 0.0) # 应用误差 theta_noisy = np.angle(Phi_diag) + phase_err beta_noisy = 1.0 + amp_err Phi_noisy = beta_noisy * np.exp(1j * theta_noisy) return np.diag(Phi_noisy) # 测试:在优化后的 Phi_opt 上叠加误差 Phi_noisy = apply_irs_errors(Phi_opt.diagonal(), phase_err_std=0.15, amp_err_std=0.03) H_noisy = H_ru @ Phi_noisy @ H_br cap_noisy = calculate_capacity(H_noisy.reshape(1,-1), W) print(f"含误差容量: {cap_noisy:.3f} bps/Hz (原始: {cap_irs_on[-1]:.3f})")提示:相位误差标准差
0.15 rad ≈ 8.6°是当前商用 PIN 二极管 IRS 的典型水平。测试表明,当phase_err_std > 0.2 rad时,容量增益下降超 40%,说明相位控制精度是 IRS 性能瓶颈,远高于幅度控制要求。
5. IRS 部署实操技巧:从仿真参数到现场调试的关键转换
5.1 仿真参数与实测设备的映射表
仿真中设置的参数必须对应真实硬件规格,否则结论无工程价值。下表列出关键参数的映射关系与典型取值:
| 仿真参数 | 物理含义 | 实测设备参考 | 典型取值 | 调试要点 |
|---|---|---|---|---|
N_irs = (8,8) | IRS 单元数 | NextNav RIS-64 | 64/128/256 | 单元数越多增益越高,但控制线复杂度平方增长 |
d = λ/2 | 单元间距 | 与频率强相关 | 28 GHz → 5.36 mm | 间距 < λ/2 可抑制栅瓣,但加工难度上升 |
bits=2 | 相位分辨率 | Metamaterial 开关 | 1-bit/2-bit/3-bit | 2-bit(4 相位)是成本与性能平衡点,3-bit 增益提升<0.5 dB |
fc=28e9 | 工作频段 | 5G n257/n261 | 26–28 GHz / 37–40 GHz | 毫米波频段路径损耗大,IRS 增益更显著;Sub-6GHz 增益有限但部署灵活 |
5.2 现场信道测量:用商用矢量网络分析仪(VNA)标定 IRS 响应
仿真中的H_br和H_ru必须通过实测校准。标准流程如下:
- 断开 IRS 控制电路,将其置于全反射状态(所有单元设为 0° 相位);
- 将 VNA 的 Port1 接基站天线,Port2 接 IRS 边缘单元馈电点(需定制探针);
- 测量 S21 参数,得到
H_br的幅度与相位(需多次测量取平均); - 同理,将 Port1 接 IRS 单元,Port2 接用户天线,测得
H_ru; - 关键步骤:更换 IRS 相位配置(如全部设为 90°),再次测量 S21,验证相位偏移是否符合预期(应接近 π/2)。
注意:VNA 测量必须在微波暗室进行,否则多径反射会污染 S 参数。若无暗室,可采用时域门(Time Domain Gating)技术滤除强反射干扰。
5.3 IRS 控制接口调试:从 MATLAB 仿真到嵌入式 MCU 的指令映射
仿真中Phi_diag是复数向量,但实际 IRS 控制器(如 STM32 或 FPGA)接收的是数字控制字。例如,某款 2-bit IRS 的每个单元由 2 根 GPIO 线控制:
| GPIO[1:0] | 相位(°) | 对应复数 |
|---|---|---|
| 00 | 0 | 1+0j |
| 01 | 90 | 0+1j |
| 10 | 180 | -1+0j |
| 11 | 270 | 0-1j |
因此,仿真输出的theta_quant需转换为控制字:
# theta_quant 单位:弧度 control_word = np.zeros(len(theta_quant), dtype=int) for i, theta in enumerate(theta_quant): if abs(theta) < 0.1: control_word[i] = 0b00 elif abs(theta - np.pi/2) < 0.1: control_word[i] = 0b01 elif abs(theta - np.pi) < 0.1: control_word[i] = 0b10 else: control_word[i] = 0b11 # 将 64 个 2-bit 字打包为 16 字节(128 bit)发送至 MCU byte_array = bytearray() for i in range(0, 64, 4): # 每 4 个单元占 1 字节 byte_val = (control_word[i] << 6) + (control_word[i+1] << 4) + \ (control_word[i+2] << 2) + control_word[i+3] byte_array.append(byte_val)现场调试时,若发现 IRS 无响应,应按顺序检查:① MCU 是否收到完整 16 字节;② GPIO 电平是否符合器件手册要求(如 1.8V vs 3.3V);③ 单元馈电网络是否虚焊。多数故障源于控制字打包错误或电压不匹配,而非算法本身。
本文还有配套的精品资源,点击获取