太赫兹UM-MIMO与IRS混合球面与平面波信道估计Matlab实践

上个月我把一套原本在毫米波频段跑得好好的信道估计代码直接搬到太赫兹参数下NMSE瞬间飙了十几个dB。起初我怀疑是信噪比设置出了问题排查到最后才发现问题出在波前模型上——当阵列孔径和通信距离进入同一量级时平面波假设已经撑不住了。今天这篇博文就围绕这个项目展开太赫兹集成UM-MIMO与IRS系统下的混合球面与平面波信道估计包括模型推导、算法选型、Matlab仿真框架以及我实际调试中踩过的一些坑。适合正在做太赫兹通信、近场信道建模或智能反射面相关课题的同行参考。1. 为什么太赫兹信道估计不能继续用平面波假设偷懒1.1 从阵列尺寸与通信距离的矛盾说起传统MIMO信道估计之所以能放心使用平面波模型是因为阵列孔径相对通信距离足够小。以28GHz、64天线为例半波长间距下阵列孔径大约0.34米瑞利距离大概是22米左右大多数室外场景都能满足远场条件波前近似成平面波带来的相位误差可以忽略。但把频率推到太赫兹情况就完全变了。以300GHz为例波长1mm64天线的线阵孔径约32mm看起来不大但瑞利距离公式是D_Ray 2L² / λ代入L0.032m、λ0.001m得到D_Ray约为2米。如果阵列为256甚至1024天线的面阵孔径做到0.1米以上瑞利距离会迅速膨胀到20米甚至更高。而太赫兹通信的典型部署场景恰恰是室内短距离、热点覆盖通信距离普遍在5到20米之间甚至更短。这意味着链路有很大概率落在瑞利距离以内波前不再是平面而是球面。这个问题的本质是频率越高波长越短同等天线数量下阵列的电气尺寸越大远场条件反而越难满足。这是一个反直觉但非常关键的结论也是整个项目选题的起点。1.2 球面波的物理本质与近场判据平面波假设下信号到达阵列各天线时波前是等相位面相位沿阵列呈线性变化。以均匀线阵为例第n个天线的接收相位可以写成φ_n -2π/λ · d · n · sinθ这里只和到达角θ有关和距离无关。而球面波模型下每个天线到辐射源的距离是真实的三维欧氏距离相位是φ_n -2π/λ · (r_n - r_ref)其中r_n sqrt((x_n - x_s)² (y_n - y_s)² (z_n - z_s)²)r_ref是参考天线到辐射源的距离。对阵列各天线的相位做泰勒展开会看到一次项对应平面波的部分二次项就是距离带来的波前弯曲也就是近场效应。远场模型实质上是把这条曲线截断成线性而近场下二次项不能忽略。判断近场还是远场业内常用两个判据一是前面说的瑞利距离二是Fresnel距离。工程上当通信距离小于D_Ray时相位误差累积到一定程度经典波束成形和信道估计都会明显退化。另外还有一点容易被忽略近场不只是一个是否问题而是一个渐变过程。距离刚好在瑞利距离附近时平面波模型误差开始显现但不是灾难性距离缩短到瑞利距离的几分之一时必须切换到球面波模型。这也引出了混合球面与平面波这个提法——在一个系统中不同链路的近场程度不同不能武断地全用球面波也不能全用平面波而应该按链路分别建模。2. UM-MIMO与IRS串联之后信道模型到底长什么样2.1 系统架构与信令链路项目采用的系统架构是基站侧部署UM-MIMO阵列例如8×8或16×16的均匀平面阵UPA天线数达到256到1024量级IRS部署在基站与用户之间常见做法是贴在墙面或楼宇外立面上包含M个无源反射单元用户端可以是单天线也可以是小规模阵列。信号路径是基站 → IRS → 用户IRS本身只有反射系数可调不具备基带信号处理能力。这正是IRS系统信道估计最难的地方我们观察到的等效级联信道是基站到IRS、IRS到用户这两段信道的乘积而不是某一段信道的直接观测。从物理传播距离来看基站到IRS的距离通常较远假设在几十米量级且基站阵列孔径虽大但相对这段距离可能还能接受这一跳倾向于平面波模型而IRS到用户的距离常常很近几米甚至一两米IRS单元数量又多例如64×64阵列孔径相对距离完全不满足远场条件这一跳必须用球面波模型。两个链路拼接在一起就构成了混合球形与平面波的物理基础。2.2 混合球面/平面波的级联信道表达为了把问题讲清楚先定义几个矩阵。设基站到IRS的链路矩阵为G维度是N×MN是基站天线数M是IRS单元数IRS到用户的信道向量为h维度M×1。IRS的反射相位配置写成一个对角矩阵Φdiag(φ)其中每个元素的幅度为1、相位可调。那么用户接收到的等效信道是h_eff h^H · Φ · G这是单用户、窄带情况下的基本表达式。实际仿真里通常假设太赫兹信道具有稀疏多径特性散射体数量有限每条路径由到达角、离开角、传播距离和复增益共同刻画。对于平面波跳段导向矢量只依赖角度比如均匀面阵的远场导向矢量 a(θ_a, θ_e) exp(-j·2π/λ·(x·sinθ_e y·cosθ_e·sinθ_a))这里的x、y是天线位置的坐标。对于球面波跳段导向矢量必须逐天线计算到用户或散射体的真实距离公式为a(r, θ_a, θ_e) [exp(-j·2π/λ·(r_1 - r_ref)), exp(-j·2π/λ·(r_2 - r_ref)), ...]其中r_i是第i个单元到目标点的距离。注意这个表达式里角度和距离耦合在一起不再像远场那样可以干净地分解。需要强调的是级联信道h_eff中IRS的相位矩阵Φ同时作用于G的列和h的行这意味着在估计过程中IRS的相位配置必须作为可设计的观测模式参与其中。每改变一次IRS的相位配置就相当于对相同的物理信道做了一次不同的投影我们正是利用多组投影来反推信道参数。2.3 从角域稀疏到角域-距离域稀疏的转变传统毫米波信道估计能够使用压缩感知根本原因是信道在角域具有稀疏性——路径数量少每条路径对应角域字典中的一个原子整个信道矩阵在角域字典下只有少量非零系数。进入近场之后球面波的导向矢量不仅由角度决定还由距离决定。同样的到达角距离不同导向矢量完全不同。因此稀疏字典必须从角域字典扩展成角域-距离域字典。这里有一个关键变化对同一个空间位置球面波导向矢量在角度维上不再是单一峰值而会在一定角度范围内出现展宽这就是近场能量扩散效应。如果继续用平面波字典去匹配会出现能量泄漏到多个角度格子上的现象稀疏性被破坏稀疏恢复算法的性能自然大幅下降。所以混合模型下的信道估计问题本质上变成了如何在一个同时包含角度和距离参数的字典上恢复稀疏系数且不同链路使用不同的波前模型。这是我决定自己写代码而不是直接套用公开OMP工具箱的核心原因——现有的大多数工具箱默认远场导向矢量改起来比重新写还要麻烦。3. 参数化信道估计的核心思路与算法选型3.1 为什么压缩感知在这里是刚需UM-MIMO阵列的天线数量动辄数百但受限于硬件成本与功耗射频链路数量通常只有个位数。也就是说我们不可能对每根天线做独立的基带观测只能通过少量射频链路和移相器网络对高维信道做低维观测。这个观测过程本身就是高度欠定的传统的最小二乘法直接失效。太赫兹信道的另一个优势是稀疏性极强。由于太赫兹频段的路径损耗严重、分子吸收显著信号能到达接收端的散射径数量非常有限室内场景下有效路径通常只有3到6条。这与压缩感知理论的要求完美匹配信号稀疏、观测欠定但可以通过稀疏恢复算法以高概率重构。然而近场带来的字典维度膨胀给压缩感知增加了难度。远场场景下二维角域字典的原子数量大约是角度网格数量的乘积例如180×180约3万个原子近场场景下还要乘上距离网格数假设距离网格取30到50个字典原子数直接到了百万量级。这个规模在Matlab里直接构造全矩阵个人电脑的16GB内存基本扛不住。因此算法设计必须兼顾恢复精度与计算复杂度。3.2 三维字典的构造与OMP恢复流程我采用的基线方法仍然是正交匹配追踪OMP但字典结构做了三维化处理。以基站到IRS链路为例假设路径数为L每条路径的未知参数包括基站侧离开角、IRS侧到达角以及对应的传播距离。字典原子的构造方式是对每个角度格点与每个距离格点的组合生成对应的球面波导向矢量再通过Kronecker积或逐列组合形成整体字典。% 生成球面波导向矢量UPA近场 function a gen_spherical_steering(lambda, pos_m, tx_pos, ref_idx) % pos_m: 阵列单元坐标维度 3 x N % tx_pos: 目标点坐标1 x 3 d vecnorm(pos_m - tx_pos); % N x 1各单元到目标点的距离 r_ref d(ref_idx); % 参考单元距离 a exp(-1j * 2 * pi / lambda * (d - r_ref)); a a / norm(a); % 归一化保持字典原子能量一致 end恢复流程分几步第一步用随机或DFT结构的组合矩阵对信道做观测得到接收向量y第二步在字典D中寻找与y相关性最强的原子第三步用最小二乘估计该原子的系数并更新残差第四步重复直到残差低于阈值或达到预设路径数最后根据恢复出的原子索引反查角度与距离参数重构完整信道。这里有一个值得注意的细节OMP字典中的相关性计算在近场字典下不能直接用二维FFT加速因为距离维破坏了傅里叶结构。如果字典规模不大直接计算内积即可如果规模很大建议把角度估计和距离估计拆成两步先用粗角度字典定位候选区域再在候选区域细化距离相当于两级搜索。3.3 降低复杂度的两个实用策略第一是分离式估计。先忽略距离维用平面波字典做一轮粗角度估计把每条路径的角度锁定到较小区间然后对每条路径在锁定角度附近生成小范围的球面波字典精细估计距离。这个策略的物理依据是近场导向矢量在角度维上虽然出现展宽但峰值位置仍然接近真实角度粗估计的偏差通常在一个波束宽度以内。实测下来两步法比直接用三维大字典的NMSE性能只差不到2dB但计算量降低了一个数量级。第二是距离网格的非均匀设计。远距离处波前曲率变化缓慢距离网格可以稀疏近距离处波前曲率变化剧烈距离网格必须加密。我采用的方案是在瑞利距离以内按平方律加密瑞利距离以外切换到线性网格并快速稀疏化。这样能把距离维的网格数量从80以上压缩到30左右而性能几乎不损失。这两个策略在Matlab项目中都是可以直接落地的后面我会给出具体的代码结构。4. Matlab仿真框架与核心代码拆解4.1 仿真参数设定与典型配置仿真参数的选择直接决定了结果是否有说服力。我从某实验室的研究工作中总结了如下典型配置大家可以根据自己的场景调整参数典型值说明载频300 GHz太赫兹代表性频段带宽2 GHz窄带假设逐子载波估计基站天线16×16 UPA256满足UM-MIMO定义IRS单元数32×321024无源反射面射频链路数4欠定观测来源路径数4角域-距离域稀疏度基站-IRS距离30 m远场链路平面波IRS-用户距离2 m近场链路球面波信噪比0~20 dB评估区间这个配置的精妙之处在于BS→IRS链路落在远场区IRS→用户链路落在近场区正好构成混合球面与平面波的典型场景。如果两个链路都改成近场模型就退化为纯球面波反而失去了混合验证的价值。4.2 球面波导向矢量生成的易错点这段代码是整个项目最容易出bug的地方我分三点说明。第一是坐标系的统一。生成导向矢量时阵列单元坐标必须与目标点坐标使用同一个世界坐标系不能一部分用局部坐标、一部分用全局坐标。我早期踩过这个坑把IRS的单元位置定义在局部坐标系里用户位置却用全局坐标导致距离计算整体偏移角度恢复了但距离参数全部错误。第二是相位参考单元的选择。球面波导向矢量有一个不可消除的公共相位选择不同的参考单元只会影响所有原子的公共相位不会影响最终恢复的信道幅度与角度参数。但如果字典构造时参考单元与信道生成时不一致会产生一个固定的相位旋转导致恢复出的信道相位失真。我的建议是统一使用阵列几何中心作为相位参考这样字典原子之间的相位关系最直观。第三是归一化。字典原子必须做能量归一化即每个导向矢量的二范数为1。否则OMP选原子时会偏向范数大的原子最小二乘的系数估计也会失真。很多初版代码NMSE差问题往往不在算法而在归一化。4.3 测量矩阵设计、字典构造与恢复代码逻辑由于射频链路数远小于天线数观测过程用一个组合矩阵W表示维度是射频链路数×天线数。实际中W由移相器网络和模拟合并器实现仿真里通常用随机相位旋转矩阵或DFT矩阵的子行。关键点在于W必须在多组IRS相位配置下重复改变才能积累足够的观测信息。% 生成测量矩阵随机相位模拟合并器 function W gen_measurement_matrix(N_RF, N, block_idx) rng(block_idx * 100 42); W exp(1j * 2 * pi * rand(N_RF, N)); W W / sqrt(N); end下面给出OMP恢复的主循环逻辑。为节省内存我采用两步法先用远场字典做角度粗估计再在锁定角度附近生成近场细字典做距离精估计。% 第一步角度粗估计远场字典 function [angle_idx, y_res] coarse_angle_omp(y, Dict_far, L) y_res y; angle_idx zeros(L, 1); for k 1:L corr Dict_far * y_res; [~, idx] max(abs(corr)); angle_idx(k) idx; % 最小二乘更新 A_sub Dict_far(:, angle_idx(1:k)); hat A_sub \ y; y_res y - A_sub * hat; end end% 第二步距离精估计近场细字典 function [r_est, gain] refine_distance(y, angle_est, tuning_args) Dict_near gen_near_field_dictionary(angle_est, tuning_args.range_grid); corr Dict_near * y; [~, idx] max(abs(corr)); r_est tuning_args.range_grid(idx); gain Dict_near(:, idx) * y; % 匹配滤波系数 end实际项目里为了获得更好的性能我在第二步之后又加了一轮联合细化把同一条路径的角度与距离参数作为初始值用非线性优化在连续参数空间里做一个局部搜索。这个步骤类似传统算法里的角度精估计但优化目标变成了球面波导向矢量与残差的内积模值。每次迭代大约几十次整体耗时可控。5. 仿真结果解读与参数敏感性踩坑记录5.1 混合模型 vs 纯平面波模型的性能对比完成仿真后我对比了三类方法一是全程假设平面波的OMP二是全程假设球面波的OMP把两条链路都建模成近场三是本文的混合模型OMP即BS→IRS用平面波、IRS→用户用球面波。在SNR15dB、路径数4、IRS-用户距离2米的场景下混合模型的NMSE约为-18dB纯平面波模型因为字典失配严重NMSE停留在-6dB左右形成明显的性能地板。这个地板正是模型误差主导的而不是噪声主导——即使把SNR推到30dB纯平面波模型的NMSE也几乎不下降这是判断模型失配的典型信号。全球面波模型在这个配置下性能接近混合模型但有一个隐患当基站到IRS的实际距离达到几十米时全球面波模型的距离网格必须覆盖很宽的范围数据量增大且距离估计的模糊问题更容易出现。混合模型的优势在于只在需要近场建模的链路上使用距离维字典降低了复杂度也减少了误判概率。5.2 我在调试中踩过的几个坑第一个坑是距离网格的尺度。最早我按线性网格从0.5米扫到50米结果近距离格点间距太大2米处的真实距离落在两个格点之间产生了不可接受的量化误差。后来改成平方律加密问题立刻消失。所以近场场景下距离网格一定要随距离非线性变化密度与波前曲率的变化率匹配。第二个坑是OMP的停止条件。路径数未知时单纯用残差阈值容易在低SNR下多恢复出虚假路径。我的经验是同时设置残差能量阈值和最大迭代次数并且在每次迭代后检查新选原子与已选原子的相关程度如果相关性超过0.9大概率是同一个物理路径的重复分裂应该终止迭代。第三个坑是IRS相位配置的随机性。最开始测试时我用完全随机的IRS相位矩阵结果发现少数几组相位下观测退化严重恢复性能波动很大。换成结构化设计——例如DFT型的相位矩阵让不同IRS单元的组合覆盖更均匀的观测空间——整体性能稳定了很多。另一个重要经验是IRS相位图案的数量要足够经验上至少要是待估计路径数的两倍以上。第四个坑与带宽有关。太赫兹频段的带宽动辄数GHz窄带假设在单载波范围内近似成立但跨子载波时导向矢量会随频率变化如果直接对宽带信号建模会引入波束偏斜误差。我在仿真中采用逐子载波处理每个子载波独立做信道估计最后在频域拼接。对于宽带系统的严格建模需要在导向矢量中加入频率依赖项那是另一个层面的话题了。最后一个坑是Matlab内存。三维字典在256天线条件下直接生成全矩阵占用内存轻松超过16GB。我的替代方案是把字典拆成按角度分段的小矩阵逐个处理或者用函数句柄代替显式矩阵在OMP内积环节逐块计算。两种方案都能把峰值内存压到4GB以内代价是循环次数增加但整体运行时间仍可接受。至于我自己的体会这个项目做完之后我对近场信道估计的认知发生了不少转变。最核心的一条是不要因为远场假设在数学上更简洁就默认它适用于所有场景。太赫兹与IRS的结合让近场不再是罕见特例而是常态。做信道估计之前先花十分钟确认每一条链路是否满足瑞利距离条件比事后调参高效得多。代码层面建议把导向矢量生成函数与恢复算法解耦这样换频率、换阵列、换波前模型时不需要动主循环只需要换一个生成函数即可。后续如果要做实时性更强的方案可以考虑把OMP换成近似消息传递或深度展开网络但那是另一个项目的故事了。