做无线物理层仿真的人,看到“太赫兹集成UM-MIMO和IRS系统的混合球面与平面波信道估计”这个标题,第一反应多半是:近场信道估计,又是个硬骨头。太赫兹频段带宽大、波长极短,配上超大规模天线和智能反射面,是6G空口Proposal期的热门组合;但当阵列口径和通信距离进入同一个数量级时,传统平面波假设会明显失配,必须引入球面波描述。这篇文章以“混合球面与平面波”为主线,把项目背后的信道模型、核心算法、Matlab实现拆开讲清楚,适合正在做太赫兹通信仿真、IRS辅助MIMO传输或者准备毕设开题的同学参考。
1. 项目背景与核心问题拆解
1.1 太赫兹频段为什么需要UM-MIMO和IRS
太赫兹频段大致覆盖0.1到10THz,是一片尚未被充分利用的频谱资源。相比5G使用的毫米波频段,太赫兹能提供几十GHz甚至上百GHz的连续带宽,理论上单链路速率可以做到Tbps级别,所以它被认为是6G超高速接入、室内短距传输、无线回传等场景最有潜力的承载频段之一。但这个频段的物理特性并不友好:自由空间路径损耗随频率平方增长,再加上水蒸气和氧气分子的吸收峰,太赫兹波能有效传播的距离通常只有几十米。想把这个频段真正用起来,必须靠大规模波束赋形把能量集中起来,UM-MIMO(超大规模MIMO)就是顺理成章的选择。
UM-MIMO的天线数通常在64到1024之间,甚至更多,可以在很小的物理尺寸上集成大量阵元。比如在1THz附近,波长只有0.3mm,半波长阵元间距约0.15mm,一片指甲盖大小的芯片上就能排布上百个阵元。这么大的阵列能形成非常窄的波束,弥补太赫兹传播损耗。但问题也随之而来:波束变窄意味着覆盖范围变小,一旦直射路径出现人体阻挡或物体遮挡,链路质量会迅速恶化。IRS(智能反射面)的价值正是在这里。它由大量低成本无源反射单元组成,不需要完整的射频链路,每个单元只控制反射相位,就能把来波重新指向用户方向,起到类似中继的补盲作用。把UM-MIMO和IRS放在一起,相当于既解决了远距离传输增益,又解决了覆盖死角,是太赫兹短距通信系统设计中很自然的技术组合。
1.2 平面波假设失效:从远场到近场
传统MIMO信道建模有一个默认前提:接收阵列离辐射源足够远,到达各阵元的波前可以近似成平面波。平面波模型下,各阵元之间的相位差只由入射角度和阵元间距决定,和辐射源到阵列的距离无关,所以阵列导向矢量可以表达成线性相位的形式,数学处理非常方便。但这个假设在太赫兹UM-MIMO场景中越来越靠不住。判断远场和近场边界最常用的标准是瑞利距离 (R_f = 2D^2/\lambda),其中 (D) 是阵列最大口径,(\lambda) 是波长。
举例来说,一个工作在1THz、包含256个阵元的均匀线阵,口径大约38mm,对应的瑞利距离约为 (2 \times (0.038)^2 / 0.0003 \approx 9.6) 米。如果阵元数提高到512,瑞利距离会膨胀到近40米。这意味着一间普通的室内会议室、一段几十米的室外回传链路,基站和用户之间的通信距离几乎都在近场范围内。近场区波前是球面波,各阵元到辐射源的距离各自不同,从不同阵元视角看过去,入射角也不再是同一个角度。此时如果还用平面波导向矢量去匹配实际信道,相位误差会随阵元序号累积,轻则降低信道估计精度,重则让波束指向完全偏掉。
1.3 混合球面与平面波模型的核心思路
如果所有链路都按球面波建模,数学和计算复杂度会大幅增加,而且实际系统里并不是所有链路都处于近场。比如基站到用户之间的直射路径可能超过几十米,对几百个阵元的小型UM-MIMO阵列来说已经进入远场区;而经IRS反射的链路距离往往很短,从用户或基站视角看IRS时仍然处于近场区。把所有路径一刀切成同一种模型,要么精度不够,要么复杂度浪费。
更合理的做法是构造混合模型:每条链路根据真实距离判断应该使用平面波还是球面波。远场链路用线性相位的平面波导向矢量,近场链路用包含距离参数的球面波导向矢量,二者共同出现在同一个信道表达式中。这个设计还有一个非常实际的优点:信道估计阶段可以把两种导向矢量拼成一个“混合字典”,稀疏恢复算法在搜索原子时会自动选出最贴合实际路径的那类原子,等价于同时完成了路径分类和参数估计。这也是项目标题里“混合”两个字落在信道估计上的直接含义。
2. 系统建模与信道数学表达
2.1 系统模型与链路关系
考虑一个典型的IRS辅助上行信道估计场景:用户终端发射导频,基站端配置 (M) 根天线的UM-MIMO阵列,IRS部署在用户与基站之间的侧面,由 (N) 个反射单元组成。每个IRS单元都通过控制器设置一个反射相位,整体形成一个对角相移矩阵 (\mathbf{\Phi} = \mathrm{diag}(e^{j\phi_1}, \dots, e^{j\phi_N}))。基站收到的信号分两路:一路是用户直接到基站的信号,用直射信道 (\mathbf{h}_d \in \mathbb{C}^{M\times 1}) 表示;另一路是用户到IRS、经IRS反射再到基站,由用户-IRS信道 (\mathbf{h}_r \in \mathbb{C}^{N\times 1}) 和IRS-基站信道 (\mathbf{G} \in \mathbb{C}^{M\times N}) 级联而成。
接收信号可以写成 ( \mathbf{y} = \sqrt{P}(\mathbf{h}_d + \mathbf{G}\mathbf{\Phi}\mathbf{h}_r)s + \mathbf{n} ),其中 (P) 是发送功率,(s) 是导频符号,(\mathbf{n}) 是噪声。注意IRS本身是无源器件,没有基带信号处理能力,它不能像有源中继那样把自己收到的信道信息解调出来。接收端能观测到的只有改变 (\mathbf{\Phi}) 后合成的级联信道 (\mathbf{G}\mathbf{\Phi}\mathbf{h}_r),所以信道估计问题的本质不是直接分别估计 (\mathbf{G}) 和 (\mathbf{h}_r),而是估计所有IRS单元相位配置下形成的一个等效信道向量。这也是IRS信道估计相比传统MIMO估计更棘手的地方。
2.2 平面波与球面波导向矢量的表达式
以均匀线阵为例,阵元间距 (d),入射角 (\theta) 定义为波到达方向与阵列法线的夹角。平面波假设下,第 (m) 个阵元相对参考阵元的相位差是 (2\pi m d \sin\theta / \lambda),导向矢量为:
[ \mathbf{a}_{\text{planar}}(\theta) = \left[1,; e^{-j\frac{2\pi}{\lambda}d\sin\theta},; \dots,; e^{-j\frac{2\pi}{\lambda}(M-1)d\sin\theta}\right]^T ]
这个式子只依赖角度,不依赖距离,所以可以预先按照角度网格构建字典,字典生成成本很低。
近场球面波模型下,每个阵元到辐射源的距离都要单独计算。假设参考阵元到源的距离是 (r),源位于与法线夹角 (\theta) 的方向上,则第 (m) 个阵元到源的距离为:
[ r_m = \sqrt{r^2 + (md)^2 - 2 r m d \sin\theta} ]
阵列响应向量中的相位项变成 (e^{-j\frac{2\pi}{\lambda}(r_m - r)}),幅度项通常可以近似忽略。很明显,球面波导向矢量是角度和距离的联合函数,不再具有 (e^{-jm\mu}) 这样的线性相位形式。因此在构造字典时,需要在角度维度之外再增加一个距离维度。距离网格一旦增加,字典规模会成倍扩大,这也是后面所有计算和内存优化的起点。
2.3 信道稀疏性与可估计性分析
在典型的太赫兹传输环境中,可分辨的多径数量远小于天线单元数。直射路径加上IRS反射路径再加少量散射簇,通常只有3到6条有效路径。这些路径在角度-距离联合域中只会占据少量原子位置,因此信道在混合字典下具有自然的稀疏性。把平面波子字典 (\mathbf{D}_p) 和球面波子字典 (\mathbf{D}_s) 拼接成总字典 (\mathbf{D} = [\mathbf{D}_p, \mathbf{D}_s]),信道向量可以写成 (\mathbf{h} = \mathbf{D}\mathbf{x}),其中 (\mathbf{x}) 是稀疏向量,非零元素的位置直接对应某条路径的到达角、距离以及该路径更适合用哪类波前模型。
稀疏性之所以重要,是因为它把信道估计问题从传统的最小二乘问题变成了稀疏恢复问题。未知参数从所有天线单元对应的信道系数,缩减成少数路径对应的复振幅和参数索引,导频开销可以从几百个符号降到几十个甚至更少。同时,可估计性还依赖于观测次数。IRS的相位配置每次改变都会产生一个新的等效信道观测,接收端搜集多次观测并将它们排列成更大的观测方程,稀疏恢复算法才能从不同相位的组合中解出支撑集。只要测量矩阵满足有限等距性质,OMP、FISTA等算法都能以很高概率恢复出稀疏信道。
3. 信道估计算法设计与实现要点
3.1 为什么传统LS/MMSE很难直接落地
如果直接把IRS到基站的信道矩阵 (\mathbf{G}) 和用户到IRS的信道向量 (\mathbf{h}_r) 全部当作未知参数估计,未知数个数会高达 (M N + N)。一个 (M=256),(N=256) 的系统,未知数就超过6.5万,LS需要导频长度不小于未知数个数,这在导频资源和时延上都不可接受。MMSE虽然可以在噪声和信道统计特性已知时获得较好性能,但它同样需要高维协方差矩阵,而且太赫兹频段的信道统计特性本身就没有成熟模型,实际系统里很难准确获取。
更重要的是,IRS信道估计中大量参与观测的是级联信道而非直射信道,直接法无法在少导频条件下把 (\mathbf{G}) 和 (\mathbf{h}_r) 分开。即便后续预编码只需要级联信道,也需要先把级联信道的稀疏表示估计出来。因此,利用信道路径稀疏性和IR前端相位的可配置性,设计压缩感知或稀疏优化算法,才是这个项目真正落地的方向。
3.2 基于OMP的稀疏信道估计流程
OMP(正交匹配追踪)是理解稀疏信道估计最直接的入口。它的核心思想是每次从字典中找到与当前残差相关性最大的一列原子,加入支撑集,然后利用最小二乘更新系数并计算新残差。算法流程如下:
- 初始化残差 (\mathbf{r} = \mathbf{y}),支撑集 (S = \emptyset),最大迭代次数 (K)。
- 计算字典所有列与残差的内积,找到绝对值最大的原子索引 (k = \arg\max |\mathbf{D}^H \mathbf{r}|)。
- 将索引 (k) 加入支撑集 (S)。
- 用支撑集对应的字典列做最小二乘更新:(\hat{\mathbf{x}}_S = \mathbf{D}_S^\dagger \mathbf{y})。
- 更新残差:(\mathbf{r} = \mathbf{y} - \mathbf{D}_S \hat{\mathbf{x}}_S)。
- 若迭代次数小于预设稀疏度且残差能量高于噪声阈值,则返回第2步,否则结束。
这段流程看起来简单,实现时有几个细节特别容易出错。字典原子必须归一化,否则幅度大的原子会在匹配步骤中占便宜,导致支撑集选错。迭代停止条件不能只用稀疏度,还要结合噪声能量,否则在高SNR下会多选出噪声原子。最小二乘更新时,如果支撑集不断增大,用 ( \mathbf{D}_S^\dagger \mathbf{y} ) 逐次求伪逆会产生大量重复计算,实测中可以用Cholesky分解增量更新,不过对入门复现可以先写简单版本。
3.3 混合模型下的原子选择与参数估计
混合字典的构造有两种常见策略。第一种是“拼接式”,把平面波字典 (\mathbf{D}_p) 和球面波字典 (\mathbf{D}_s) 直接横向拼接,交给OMP统一搜索。优点是一个算法框架同时解决两种波前模型,不需要提前判断每条路径属于远场还是近场;缺点是字典规模膨胀,内存和搜索开销大。第二种是“分级式”,先用平面波字典粗估计,找出能量不够集中或残差偏大的区域,再在这些区域引入球面波距离网格做细估计。该方法在理想情况下计算效率更高,但逻辑分支较多,对初学复现并不友好。
我在实际项目中的建议是:阵列规模不大、IRS单元数在256以内时,直接拼接字典最省事,多花一点内存换取程序清晰度和算法稳健性;当天线数上千后,再考虑分块字典或树形搜索。有一点要特别注意,球面波子字典的距离网格不能等间隔取到很密,否则相邻距离的原子相似度太高,OMP容易出现原子混叠。通常距离分辨率选0.5到1米即可满足大多数室内近场场景。
3.4 性能评估指标与对比维度
评价信道估计算法好坏,不能只看一条NMSE曲线。最常用的指标是归一化均方误差:
[ \mathrm{NMSE} = \mathbb{E}\left[\frac{|\hat{\mathbf{h}} - \mathbf{h}|^2}{|\mathbf{h}|^2}\right] ]
用来衡量估计信道和真实信道之间的差异。但通信系统最终关心的是可达速率,仅看NMSE可能误导,因为某些对波束成形影响不大的分量误差不会显著改变速率。更完整的对比应当同时给出“NMSE对SNR曲线”“NMSE对导频数曲线”“可达速率对SNR曲线”三组结果。此外,还要记录算法运行时间和峰值内存,尤其当混合字典非常大时,复杂度是否能被接受往往决定方案能否实用。
| 指标 | 计算方式 | 关注点 |
|---|---|---|
| NMSE | 估计值与真实值差的能量比 | 信道估计基本精度 |
| 可达速率 | 基于估计信道做波束成形后的频谱效率 | 系统最终性能 |
| 支撑集正确率 | 恢复出的稀疏位置与真实路径重叠比例 | 稀疏重构有效性 |
| 运行时间 | 单次估计耗时的平均值 | 算法复杂度 |
| 峰值内存 | 字典矩阵和相关变量的占用 | 工程实用性 |
4. Matlab源码结构与实操要点
4.1 源码模块划分
我的习惯是先把一个Matlab项目拆成可独立调试的模块,而不是把所有代码堆在一个main脚本里。以这个信道估计项目为例,至少需要以下功能模块:
main.m % 主脚本:参数设置、蒙特卡洛循环、画图 setupParameters.m % 载波频率、阵列规模、IRS单元数、SNR等 generateChannel.m % 依据混合波前模型生成真实信道 buildDictionary.m % 构造平面波+球面波混合字典 estimateChannel.m % OMP或改进稀疏恢复算法 plotResults.m % NMSE、可达速率曲线绘制这种拆分的好处是改参数、换算法都不会牵动底层代码。值得注意的问题是变量维度必须统一。我见过很多复现代码报错,都是因为某个函数内部把 (M \times N) 的矩阵转置后传给了下一个函数,维度没对齐。建议在函数开头用注释写明输入的维度,比如% G: M x N、% h_r: N x 1,调试时会节省大量时间。
4.2 关键代码片段剖析
字典构建是最核心、也最容易写错的函数。平面波字典只需要遍历角度网格,球面波字典需要在角度网格基础上嵌套距离网格:
function D = buildDictionary(lambda, d, M, angles, distances) Dp = zeros(M, length(angles)); for i = 1:length(angles) Dp(:, i) = planarSteering(lambda, d, M, angles(i)); end Ds = zeros(M, length(angles) * length(distances)); col = 0; for i = 1:length(angles) for j = 1:length(distances) col = col + 1; Ds(:, col) = sphericalSteering(lambda, d, M, distances(j), angles(i)); end end D = [Dp, Ds]; end球形导向矢量函数里,核心是求解每个阵元到源的距离。推荐先按公式循环写,验证正确后再用向量化方法替换。循环版本的代码可读性强,不容易出维度错误。OMP估计函数同样可以用一个简洁的循环实现。
function x_hat = omp(D, y, K, noise_thresh) r = y; S = []; x_hat = zeros(size(D, 2), 1); for iter = 1:K corr = D' * r; [~, idx] = max(abs(corr)); S = union(S, idx); x_hat(S) = D(:, S) \ y; r = y - D(:, S) * x_hat(S); if norm(r)^2 < noise_thresh break; end end end这个函数每次迭代都用D(:, S) \ y重复做最小二乘,效率不算高但够清晰。要加速的话,可以在支撑集扩张时利用上一次的中间结果做增量更新,但初学时建议先保证正确。还有一点,如果字典是复数矩阵,D'是共轭转置,不能写成D.',这两种转置在复数内积上有本质区别,非常容易埋坑。
4.3 运行环境与依赖
Matlab版本建议使用R2020a以上,核心代码全部用矩阵运算实现,不依赖特定工具箱也能运行。生成调制信号、计算误码率时会用到Communications Toolbox,但信道估计主体部分用标准Matlab即可。如果机器内存只有8GB,需要严格控制字典规模。一个 (M=256),角度网格720个,距离网格100个的系统,球面字典列数有7.2万,加上平面字典,按复数double计算会占用近300MB内存,看似不大,但后续相关运算中临时矩阵会成倍放大内存占用。这时候先用128天线、64个IRS单元跑通,再逐步放大规模是最稳妥的路线。
4.4 仿真参数配置建议
下面这组参数适合作为复现项目的基线配置,数据量适中,能较快验证算法趋势:
| 参数 | 推荐取值 | 说明 |
|---|---|---|
| 载波频率 | 1 THz | 波长0.3mm,近场效应明显 |
| 基站天线数 | 128 | 入门级UM-MIMO规模 |
| IRS单元数 | 64 | 级联信道规模可控 |
| 阵元间距 | 0.5λ | 均匀线阵常用设置 |
| 距离范围 | 2米到20米 | 覆盖近场区域 |
| 角度网格 | 1° | 360个角度原子 |
| 距离网格 | 0.5米 | 避免字典过大 |
| 多径数 | 3 | 保证稀疏性 |
| 信噪比 | 0~20 dB | 观察误差变化 |
| 蒙特卡洛次数 | 100 | 足够平滑曲线 |
实际调试时,先固定信噪比为20dB,把随机种子固定,逐次打印估计出的角度、距离和复振幅,确认单条路径已经恢复正确,再放开参数做完整模拟。盲改参数会让自己失去判断依据。
5. 常见问题与排查技巧实录
5.1 内存不足与矩阵维度爆炸
这是复现近场信道估计最常遇到的坑。混合字典动辄几万列,再加上复数运算,内存会非常紧张。我见过有人把距离网格设成0.1米,一下生成几十万个原子,程序直接卡死。解决思路有几个方向:把字典按距离切块,每次只加载一部分原子做匹配,选择后再合并支撑集;或者用single代替double存储字典,内存减半但对性能影响不大;也可以利用均匀线阵的Vandermonde结构,用FFT快速计算导向矢量相关值,避免显式生成完整字典。综合来说,先粗后细的两级字典是内存和精度之间最好的折中。
5.2 估计误差居高不下的排查思路
如果NMSE曲线很平或下降缓慢,先不要怀疑OMP实现,优先检查信道生成与字典是否匹配。常见问题有三个:球面波导向矢量中的参考距离用错,导致相位项出现系统性误差;用户距离恰好落在近场/远场边界附近,算法用平面波原子去拟合球面波路径,能量泄漏到多个相邻原子;IRS相位配置存在量化误差,实际观测矩阵与理论矩阵不一致。定位方法很简单,固定一条直射路径,SNR设为30dB,打印估计出的原子索引与真实参数,如果索引偏了一格,就回到字典构造函数中逐项核对数组顺序。
5.3 近场/远场切换边界的选择错误
混合模型中,到底用平面波还是球面波,不能完全依赖瑞利距离公式做硬切换,否则边界附近会出现明显的性能跳变。我的做法是设置一个过渡带:距离小于0.3倍瑞利距离时强制用球面波原子;距离大于3倍瑞利距离时强制用平面波原子;中间区域同时保留两种原子,让OMP通过内积最大原则自行选择。这样虽然会增加一部分字典列数,但换来了边界处的稳定性能,很值得。还要注意瑞利距离计算用的是阵列最大口径,不是天线阵元总数,计算时要先换算成物理尺寸。
5.4 运行时间过长与算法加速
OMP在大字典上的计算瓶颈在每轮做全字典内积。可以尝试用Kronecker积分解二维导向矢量,把角度和距离搜索拆成两个一维搜索,大幅减少搜索量;改用FISTA或ADMM这类阈值迭代算法,避免逐次选原子的顺序操作;或者在SNR点之间用parfor并行跑蒙特卡洛,提高整体吞吐。对课程设计或毕设来说,没必要追求极致速度,先把曲线趋势跑出来,再讨论优化。
6. 个人实践经验与后续扩展
最后聊一点我实际做近场信道估计项目的体会。不要一上来就搭“全球面波”模型,也不要用纯平面波假装近场不存在。最稳妥的路线是先用中等规模的阵列把近场/远场混合情况跑通,观察哪些路径在平面波模型下误差大,再用混合字典把这类路径“救回来”。复现这类源码时,我强烈建议固定随机种子,把信道生成和字典构建单独拉出来调试,确认同一组参数下字典原子与信道相位完全吻合后,再进入估计算法部分。否则一旦全流程出错,很难分清是信道模型的问题还是算法实现的问题。
后续扩展方向其实很明确:把IRS相移的码本量化误差纳入观测模型,让算法直接估计非理想硬件下的级联信道;或者用学习式稀疏重构网络替代OMP,在保持稀疏性的同时进一步降低导频开销;再往上还可以结合波束对齐和信道估计做联合设计,减少训练阶段的开销。做通信仿真有个绕不开的常识:先让理想模型得到准确解,再谈非理想因素。这个原则在太赫兹混合信道估计项目中非常重要。