太赫兹UM-MIMO与IRS混合信道估计：球面波与平面波联合稀疏恢复

最近这个太赫兹集成UM-MIMO和IRS系统的混合信道估计项目在仿真圈子里讨论度挺高，版本编号都出到14942期了。我也照着思路自己完整跑了一遍，把代码结构、信道建模、字典设计这些核心环节都重新捋清楚了。这个项目本质上不是单纯调一个函数就能出结果的demo，它背后牵扯到太赫兹频段的路径损耗处理、超大规模天线阵列的球面波假设、智能反射面的相位调控，以及一套能同时处理近场和远场分量的稀疏估计算法。如果你正在做6G相关的仿真预研，或者准备把IRS和超大规模MIMO结合起来做信道估计，这篇内容应该能给你省不少时间。

先简单交代一下这套系统要解决什么问题。太赫兹频段频率高、波长短，天线尺寸可以做得很小，所以能在有限面积里塞下成百上千个阵元，这就是UM-MIMO（超大规模MIMO）。天线多了之后，目标距离相对于阵列孔径就不再能看作无穷远，传统平面波假设开始失效，必须考虑球面波。同时，端到端直射路径在太赫兹频段通常很弱，中间得放IRS做反射辅助通信，反射链路又进一步引入了近场特性。这样一来，接收信号里同时混合着平面波分量和球面波分量，估计难度直接上了一个台阶。

1. 项目概述：这套系统到底在估计什么信道

1.1 太赫兹+UM-MIMO+IRS：为什么偏偏是这三个词

太赫兹通信的定位是6G时代的核心候选技术，频段大概在0.1 THz到10 THz之间，可用带宽大、速率高，但代价是路径损耗非常严重。高频带来的另一个问题是波长在毫米级别甚至更小，同样的物理尺寸可以排布更多天线单元，所以UM-MIMO和太赫兹几乎天然绑定。天线数量一旦到几百上千，阵列孔径可能达到几十个波长，这时候“目标在远场”这个假设就站不住了，必须区分近场和远场。

IRS则是用来改善覆盖的，它由大量无源反射单元组成，每个单元可以独立调节相位，把入射信号朝指定方向反射。把IRS放进太赫兹系统里，能显著提升反射链路增益，但它也带来了一个新的麻烦：反射路径的几何关系更复杂，信道参数里同时包含发射端到IRS、IRS到接收端这两段子信道，每一段都可能是近场或远场，需要联合估计。

1.2 “混合球面与平面波”到底在说什么

我一开始看这个项目标题时，也以为“混合球面与平面波”只是某个算法名称，后来推了一遍模型才发现，这是对信道物理特性的一个很准确的描述。

近场条件下，阵列接收到的信号波前是曲面，不同阵元到同一个点的距离差不再是一个简单的线性函数，而是二次函数甚至更高阶，所以导向矢量里会引入距离相关的相位项。远场条件下，波前近似成平面，导向矢量只跟到达角度或离开角度有关。太赫兹系统里，天线数量大、阵列孔径大，而通信距离可能在几米到几十米之间，所以有些链路处于近场，有些链路处于远场，甚至同一条链路里，部分阵列单元对目标点处于近场、另一部分处于远场。因此，信道建模时不能一刀切，得用一个能兼容两种情况的混合模型，把每条路径划分为多个子路径，子路径再根据距离和瑞利距离的关系选择球面波或者平面波表达式。

这套说法听起来抽象，但实际仿真时就是判断一个瑞利距离：

R_D = 2D^2 / λ

其中D是阵列孔径，λ是波长。如果通信距离小于R_D，用球面波模型；远大于R_D，就用平面波模型。天线规模大了之后，R_D可能达到几十米，所以通信距离处于过渡区的情况非常常见，这就是为什么要做“混合”模型。

2. 系统建模与信道模型设计

2.1 几何结构与参数设定

搭建这个仿真，第一步是定几何参数。我参考项目里的常见做法，设置发射端为UM-MIMO均匀面阵，阵元数是64×64，也就是4096个天线；接收端是32×32的均匀面阵；IRS是64×64的反射单元。

工作频率取300 GHz，波长就是1 mm。如果阵列边长是64个单位波长，那么阵列孔径约为0.064 m，瑞利距离大约为2×(0.064)^2 / 0.001 ≈ 8.19 m。我在仿真里让发射端到IRS的距离是5 m，IRS到接收端是4 m，都小于瑞利距离，所以这两段子信道都要按近场球面波处理。直射路径是12 m，大于瑞利距离，按远场平面波处理。这样就构成了一个自然的混合场景。

需要说明的是，这里距离设置不是随便拍的，是为了让近场和远场分量同时出现。如果你把距离全部调到30 m开外，整个模型退化成纯平面波，就失去混合估计的意义了。

表1是我用的一组基准参数。

2.2 球面波和平面波的导向矢量怎么推

这是整个信道路径损耗之外最核心的数学环节。平面波导向矢量很简单，以均匀线阵为例，第n个阵元的相位是：

a_n(θ) = exp(-j 2π n d sinθ / λ)

其中d是阵元间距，θ是到达角或离开角。均匀面阵就是两个维度分别算之后做Kronecker积。这个大家应该都熟。

球面波导向矢量要麻烦一点。假设目标点坐标为p = [x0, y0, z0]，阵列第m个阵元的坐标为p_m = [x_m, y_m, z_m]，那么目标到第m个阵元的距离是：

r_m = ||p - p_m||_2

球面波导向矢量的第m项是：

g_m = (r_ref / r_m) * exp(-j 2π r_m / λ)

这里r_ref是参考距离，通常取阵列中心到目标的距离，用来归一化幅度。对比一下就能发现，平面波只考虑角度，球面波把每个阵元到目标点的真实距离差异全部算进去了，所以它能精确刻画近场的距离-角度耦合。

关键点在于，当距离r_m很大时，r_m对角度求近似可以展开成r_ref + d sinθ + 二阶项，二阶项可忽略时球面波自动退化成平面波。所以“混合模型”不是两套独立表达式硬拼，而是用球面波作为通用形式，在远场条件下自然退化。不过实际实现时，为了降低字典规模，还是会分别生成近场子字典和远场子字典，再拼接起来。

2.3 接收信号模型与导频设计

系统的接收信号可以写成：

y = H_total * s + n

其中s是发送导频符号，n是加性高斯白噪声。这里的H_total是级联信道，由直射信道和IRS反射级联信道组成。

IRS反射信道有一个非常重要的结构特性，就是它通常是三个矩阵的乘积：

H_irs = H_R * Φ * H_T

H_T是发射端到IRS的信道矩阵，H_R是IRS到接收端的信道矩阵，Φ是一个对角矩阵，对角元素是IRS每个单元的反射系数。反射系数的幅度通常设为1，相位由IRS的调控算法决定。

导频设计方面，我在仿真里给每个发射天线分配了正交导频序列。UM-MIMO天线数量太多，完全正交导频的时频开销是天线数目级别的，实际不可行。所以项目里通常采用压缩感知的思路，让多个天线共享少量导频，然后利用信道在角度域和距离域的稀疏性恢复完整信道。这也是后面为什么用OMP等稀疏恢复算法的主要原因。

3. 信道估计算法核心思路与MATLAB实现

3.1 为什么不上LS，而是用稀疏恢复

最开始我也想过，直接最小二乘估出来不就完了。但算一笔账就明白了：发射天线N_t是4096，接收天线N_r是1024，信道矩阵H的维度是1024×4096，未知参数超过四百万。而导频符号如果只有16个，接收信号总维度是接收天线数×导频数，也就是1024×16，只有一万六千多个观测值。这是一个严重欠定问题，LS解不唯一，噪声稍微大一点估计就崩了。

所以必须利用信道的稀疏性。太赫兹信道多径数量很少，工程上一般认为散射路径在3到10条之间，每条路径对应一个角度和距离参数。也就是说，信道矩阵在大字典下只有少量非零系数，这就天然满足稀疏恢复的前提。用OMP这类贪婪算法加上字典矩阵，就能从少量观测里恢复出高维信道。

3.2 字典是怎么生成的

字典是整个算法的灵魂。我这里的做法是生成一个由原子组成的过完备矩阵，每个原子对应一组参数组合，包括发射角、到达角、距离参数。字典的维度是观测维度乘以原子数量，原子数量由网格密度决定。

下面是MATLAB生成球面波原子的核心代码，我做了一些注释。

function dict = gen_spherical_dict(x_grid, y_grid, z_grid, array_pos, lambda, ref_dist) % 输入: % x_grid, y_grid, z_grid: 空间采样网格 % array_pos: 阵列阵元坐标, 维度为 N x 3 % lambda: 波长 % ref_dist: 参考距离 % 输出: % dict: 各候选位置对应的导向矢量集合 N = size(array_pos, 1); num_grid = numel(x_grid); dict = zeros(N, num_grid); for g = 1:num_grid p = [x_grid(g), y_grid(g), z_grid(g)]; r = sqrt(sum((array_pos - p).^2, 2)); dict(:, g) = (ref_dist ./ r) .* exp(-1j * 2 * pi * r / lambda); end end

注意这里用了三层循环的思维展开，实际运行会非常慢。我自己跑的时候做了向量化处理，把坐标网格转成矩阵一次性计算距离，速度能提升几十倍。建议你自己写代码时也优先做矩阵运算，不要逐网格循环。

平面波字典更简单，本质就是常规DFT矩阵的一部分：

function dict = gen_planar_dict(theta_grid, array_pos, lambda) % 只考虑角度, 不考虑距离 N = size(array_pos, 1); num_angles = length(theta_grid); dict = zeros(N, num_angles); for a = 1:num_angles dict(:, a) = exp(-1j * 2 * pi * array_pos(:,1) * sin(theta_grid(a)) / lambda); end end

3.3 OMP迭代实现细节

OMP算法大家应该都写过，这里只说在IRS系统里特别要注意的两个点。

第一，级联信道的字典不是发射字典和接收字典的简单拼接，而是它们的Kronecker积。具体来说，IRS反射链路里，每个路径参数对应一个从发射子字典选的原子和一个从接收子字典选的原子，联合原子是这两个导向矢量的外积的向量化。这个外积计算量很大，我建议预先计算每个单独原子的内积表，避免在迭代里重复计算。

第二，OMP的残差更新要作用在整个观测向量上，不是逐块更新。有些同学把观测矩阵拆成子块分别做OMP，结果恢复出来的信道拼接后不一致，NMSE很差。原因是级联信道有跨路径的结构相关性，分开估计丢失了这种相关性。

下面是一个OMP主循环的骨架。

r = y; % 残差 support = []; for iter = 1:max_iter proj = A' * r; % A是字典矩阵 [~, idx] = max(abs(proj)); support = union(support, idx); % 最小二乘更新系数 x_est = A(:, support) \ y; r = y - A(:, support) * x_est; if norm(r) < tol break; end end

这里的A是字典矩阵，维度是观测数×原子数。实际操作中A通常没办法完整存下来，我会用函数句柄代替显式矩阵，让A的乘法通过快速变换实现。这种“隐式字典”能大幅减少内存占用，后面详细说。

4. 仿真结果与性能分析

4.1 NMSE随信噪比的变化

我跑完基准参数后，统计了信道估计的归一化均方误差（NMSE），定义是：

NMSE = ||H_est - H_true||_F^2 / ||H_true||_F^2

结果是在0 dB时NMSE大约在-5 dB左右，15 dB时降到-25 dB左右，20 dB以上能到-30 dB以下。这个趋势符合预期，说明稀疏恢复算法在高信噪比下能逼近理论的Cramer-Rao界。

值得说明的是，这里直接看NMSE绝对值意义不大，更关键的是对比“纯平面波字典”和“混合字典”的差异。我在同样条件下用纯平面波字典做估计，15 dB时的NMSE比混合字典差了将近10 dB，原因很简单——近场路径被平面波字典无法表示，能量泄漏到其他原子上了。

4.2 IRS单元数对性能的影响

IRS单元数从16×16增加到64×64，性能提升非常明显。主要原因是IRS单元越多，反射链路增益越强，等效信噪比越高。但要注意，单元数增加会导致级联信道的字典维度爆炸式增长，仿真时间可能从几分钟涨到几十分钟。我实际测试时发现，32×32的IRS单元配合256×256的发射天线，单次蒙特卡洛仿真CPU时间大约在40秒左右；升到64×64后，单次要超过5分钟。如果你只是验证算法流程，建议先从32×32开始跑通再上规模。

表2是一次典型仿真的运行时间记录。

4.3 复杂度分析

复杂度主要来自字典矩阵乘法和OMP迭代中的最小二乘求解。假设原子数为P，观测维度为M，OMP迭代K次，那么每次迭代中A'*r的计算复杂度是O(KMP)。在P比较大时，这部分是绝对瓶颈。

我做过一个分析：当P=50000，M=4000，K=8时，A'*r每次大约要执行2×10^8次复数乘法。在普通笔记本上大约需要0.2到0.5秒。如果蒙特卡洛次数是1000，光这一步就要好几分钟，所以实际仿真一定要做降采样或者并行化。

5. 常见问题与排查技巧

5.1 内存爆炸：字典矩阵存不下来

这是做UM-MIMO仿真最容易踩的坑。假设观测维度M=4096，原子数P=50000，如果直接构造满秩复数矩阵A，需要存储的数据量是4096×50000×16字节，大约3.3 GB。这还只是单精度，如果默认double，直接翻倍到6.6 GB。我一开始没注意，直接在服务器上跑，结果OOM了。

解决办法有两个。第一个是通过函数句柄实现隐式字典，不显式存储A。第二个是把矩阵分块，比如每次只处理部分原子列，算出内积之后累加。我后来采用的是隐式字典方案，在内存占用上从6 GB降到了几百 MB。

还有一个很实用的技巧：生成字典时先把所有向量归一化。归一化可以提前算好，不要每次调用A'*r都重新算模长，无谓增加计算量。

5.2 OMP迭代不收敛或者收敛慢

最典型的原因是网格分辨率太低，真实路径角度落在两个网格点之间，能量被分散到相邻几个原子，导致OMP一次选不准。我的经验是先用粗网格快速估计路径数，再在估计值附近做一次细网格局部搜索，性能和速度都能兼顾。

另一个原因是残差阈值设置得太严格。如果噪声比较大，残差下降到一定程度后基本就是噪声地板，继续迭代不仅不提升性能，还会把噪声拟合成信号。我通常的阈值设成接收信号总能量的1%到5%，具体根据SNR调整。

5.3 MATLAB版本和工具箱问题

这个项目用到的函数基本都是MATLAB基础语法，不需要额外的通信工具箱。但我建议用2021a之后的版本，有几个原因：新版对复数矩阵乘法的优化更好，且函数句柄和稀疏矩阵处理速度有提升。如果你用的是旧版本，运行时间可能会差两倍以上。

另外，如果要用parfor做蒙特卡洛并行，注意随机数流的控制。每轮仿真之前先重置随机数生成器，否则并行池里每个worker有可能生成相同的噪声序列，结果会完全失真。我以前就因为没设置随机数种子，跑了1000次蒙特卡洛，结果所有样本都一样，查了半天才发现问题。

6. 一点个人经验

这个项目我前前后后调了两周，最耗时间的不是写代码，而是理解球面波和平面波在这个系统里怎么平滑过渡。你如果也是刚开始接触，我建议先别急着上大阵列，把16×16的天线和16×16的IRS跑通，然后用一个已知的入射方向去验证字典和恢复算法是否准确，确认无误再放大规模。

另外强烈建议把中间变量全部保存下来，尤其是字典矩阵和真实信道，这样每次改动参数后可以直接对比，不用从头跑仿真。我这边最后整理了一套可复用的流程：先生成几何参数，再生成字典，然后做OMP估计，最后画NMSE曲线。整套逻辑理顺之后，换参数场景基本只需要改开头几行配置。

如果你想在这个基础上继续扩展，可以考虑把IRS的相位优化和信道估计联合起来做，或者用深度学习网络替代OMP做近远场混合估计。现在这个框架的好处是模块化清晰，替换算法时不需要动信道生成部分，这也是我做这个项目最大的收获。