做无线通信仿真的人都知道信道建模要是“不够像回事”后面算出来的性能指标基本就是自嗨。最近我在做一版无人机通信系统仿真核心是把空中智能反射面IRS挂到无人机群上并且把信道从常见的瑞利/莱斯升级成 Nakagami-m 多径衰落 逆伽马阴影衰落 的复合模型。这个“无人机群改型”的项目用 MATLAB 从零写了一套蒙特卡洛链路仿真最终能输出中断概率、遍历容量、误码率等关键指标。整个过程踩了大大小小不少坑尤其是逆伽马阴影分布怎么生成、IRS 的相移怎么和无人机群的拓扑配合网上能直接照搬的代码很少。这篇文章就完整复盘一遍模型怎么设计、公式怎么落地、MATLAB 怎么实现、结果怎么看以及哪些细节最容易栽跟头。1. 项目概述与核心动机1.1 从单无人机到无人机群这个“改型”到底改了什么最初这版代码处理的场景很朴素一个基站、一个无人机、无人机上挂一块智能反射面给地面盲区用户做辅助通信。这种单无人机 IRS 的好处是部署灵活高度上有天然优势比地面 IRS 更容易避开障碍物。但真正跑仿真的时候会发现单块 IRS 的覆盖范围有限反射链路一旦遇到强阴影衰落性能掉得特别快。所以我把模型扩成了“无人机群 多块空中 IRS”多架无人机各自携带反射面在目标区域上方形成分布式反射阵列。用户同时收到来自多个无人机 IRS 的反射信号接收端可以等效成“多个空中反射面协作”。这里的核心变化不只是把无人机数量从 1 改成 N还引入了协同问题。每块空中 IRS 都有自己的位置、高度、与用户之间的信道状态反射信号的相位对齐和功率叠加方式和单无人机完全不同。在 MATLAB 里这意味着信道矩阵的维度从单个 IRS 的 M 个反射单元变成了 K 架无人机、每架 M 个单元的 K×M 结构而且每个单元对应的路径损耗、小尺度衰落、阴影衰落都要单独生成。整个“改型”相当于把单链路仿真升级成分布式协同仿真这比简单重复 K 次单无人机仿真要复杂得多因为合并信号的统计特性变了。1.2 为什么偏偏选 Nakagami-m 和逆伽马阴影很多教材级仿真默认用瑞利信道因为它解析简单、数学形式漂亮。但实际空对地链路里直视径可能出现也可能被遮挡反射的散射体分布很不均匀用瑞利一个参数根本描述不过来。Nakagami-m 分布多了一个形状参数 m能很好地统一“轻度衰落”和“重度衰落”两种情况m1 时退化成瑞利m 越大越接近高斯/视距信道。对无人机高度变化、飞行环境从城区到郊区切换的场景来说用一个 m 值就能调整信道恶劣程度非常实用。真正让我觉得这个项目有价值的是同时引入了逆伽马分布的阴影衰落。传统阴影衰落常用对数正态分布但它的尾巴偏“短”对极端恶劣阴影的描述不够充分。逆伽马分布有更明显的重尾特性当时在文献里看到用逆伽马建模阴影衰落时我第一反应是这不就是“极端环境友好型”阴影模型吗它能模拟某一瞬间反射路径被严重遮挡、信道功率骤降的极端情况。把 Nakagami-m 的小尺度衰落和逆伽马的大尺度阴影乘在一起得到的复合信道比单层衰落模型要贴近无人机飞行时的真实信道变化这也是这个项目区别于其他普通无人机通信仿真的核心卖点。2. 信道建模与系统设计2.1 系统拓扑与信号模型整个系统我画了一个最简拓扑地面基站 TX 发射信号经过无人机群上的 K 块空中 IRS 反射最终到达地面用户 RX。每块 IRS 上有 M 个反射单元基站到每架无人机、每架无人机到用户的下行链路由直视路径和散射路径叠加而成。信号经过每一块空中 IRS 后反射单元会对入射信号附加一个可调的相移。理想情况下所有反射单元在用户端都做到同相叠加用户收到的等效信号就是 K×M 个反射分量的相干求和。用公式表示用户端的基带等效接收信号可以写成y sqrt(P_t) * h_eq * x n其中 P_t 是发送功率x 是发送符号n 是零均值复高斯噪声h_eq 是等效信道系数。在没有 IRS 相位优化时h_eq 是每个反射链路的随机叠加在理想相位优化后等效信道幅度变成所有反射分量幅度的代数和|h_eq| Σ_{k1}^{K} Σ_{m1}^{M} |h_{k,m}|这里的 |h_{k,m}| 是第 k 架无人机第 m 个反射单元的等效信道幅度。如果进一步考虑路径损耗、阴影衰落和小尺度衰落每个分量可以建模为三个因子的乘积路径损耗幅度因子、阴影衰落幅度因子、Nakagami-m 小尺度衰落幅度因子。用数学表达式写就是|h_{k,m}| sqrt(L_{k,m}) * sqrt(S_{k,m}) * |g_{k,m}|其中 L_{k,m} 是路径损耗功率因子S_{k,m} 是逆伽马阴影功率因子|g_{k,m}| 是 Nakagami-m 小尺度幅度。合成后的信噪比就是γ P_t * |h_eq|² / N_0这个模型的好处是把“无人机群改型”的关键点都包含进去了无人机数量 K 决定反射支路数IRS 阵元数 M 决定每块反射面的增益Nakagami-m 参数 m 决定小尺度衰落的剧烈程度逆伽马参数决定阴影有多“毒”。2.2 Nakagami-m 衰落不只是“广义瑞利”Nakagami-m 分布的概率密度函数长这样f_g(r) (2 m^m r^{2m-1}) / (Γ(m) Ω^m) * exp(-m r² / Ω)其中 Ω 是平均功率也就是 E[|g|²]m 是衰落严重程度参数。m 越小衰落越重m0.5 时比瑞利还差m1 就是瑞利m→∞ 趋于无衰落的静态信道。在无人机空对地场景中低空飞行且遮挡多时 m 取 0.8~1.2高空视距条件好时可取 3~5。MATLAB 里生成 Nakagami-m 幅度有一个非常方便的关系如果功率增益 W |g|² 服从 Gamma(m, Ω/m) 分布那么 |g| sqrt(W) 就服从 Nakagami-m 分布。Gamma(m, Ω/m) 的均值为 Ω方差是 Ω²/m。参数 m 越大W 的方差越小意味着信道越稳定。因此在仿真里我只需要用gamrnd(m, Omega/m)生成功率增益再开根号就是幅度W gamrnd(m, Omega/m, K, M); % 每个IRS单元的小尺度功率增益 g_amp sqrt(W); % 每个单元的信道幅度服从Nakagami-m这个生成方式比直接逆变换采样要快得多而且天然支持任意正实数 m不需要做特殊函数求逆。2.3 逆伽马阴影被忽略的重尾效应阴影衰落用来描述大尺度遮挡导致的中值功率波动传统上常用对数正态分布。但逆伽马分布有一个对数正态没有的优点解析上方便和 Gamma 分布做乘积而且在参数选择合适时重尾更明显能模拟“偶尔一下信道特别差”的突发事件。对无人机群协同来说风速变化、机身姿态变化、障碍物遮挡切换都可能造成这种突发性阴影所以逆伽马是有实际意义的。逆伽马分布的概率密度函数定义为f_S(s) β^α / Γ(α) * s^(-α-1) * exp(-β/s)其中 α 是形状参数β 是速度参数也常被当作尺度参数处理。这个分布和 Gamma 分布互为倒数关系如果 Y 服从 Gamma(α, rateβ)那么 S 1/Y 服从逆伽马(α, β)。这里的“速率参数”要注意和 MATLAB 的gamrnd的尺度参数区分开gamrnd(a,b)返回 Gamma(a, b)其中 b 是尺度参数所以我们需要用gamrnd(a, 1/b)来模拟速率为 b 的 Gamma 分布。在仿真里生成逆伽马阴影就是这样一行代码S 1 ./ gamrnd(alpha, 1/beta, K_links, 1); % 生成逆伽马阴影功率因子为什么说它重尾因为逆伽马分布的尾巴按 s^{-α-1} 衰减当 α 比较小比如 2~4时尾巴衰减得很慢仿真中常会冒出一个特别大的 S导致某条反射链路瞬时信噪比异常低。这其实很符合低空无人机场景一阵风把机身吹偏或者一栋楼突然插到反射路径里信道就会瞬间恶化。3. MATLAB 仿真实现从链路级到系统级3.1 仿真主脚本与参数体系这个项目我采用蒙特卡洛仿真框架不推导复杂的闭合表达式而是直接生成大量随机信道样本对每个样本计算瞬时信噪比再统计低于门限的比例作为中断概率。主脚本的第一步是把所有系统参数模块化方便单独研究某个参数的影响。参数体系分成三块通信基本参数、无人机与 IRS 配置参数、信道统计参数。我在代码里是这样初始化这些参数的%% 通信基本参数 Ptx_dBm 23; % 发射功率 23 dBm N0_dBm -104; % 噪声功率 -104 dBm bandwidth 1e6; % 带宽 1 MHz %% 无人机与IRS配置 K 4; % 无人机数量 M 16; % 每块IRS反射单元数量 height 100; % 无人机高度米 dist_BS_UAV 300; % 基站到无人机距离米 dist_UAV_UE 500; % 无人机到用户距离米 %% 信道统计参数 PL_dB -80; % 标称路径损耗包含距离衰减 m 2; % Nakagami-m 形状参数 Omega 1; % Nakagami-m 平均功率 alpha 4; % 逆伽马形状参数 beta 3; % 逆伽马速度参数 gamma_th_dB 0; % 中断门限dB Nmc 1e5; % 蒙特卡洛次数这里有个经验路径损耗我一开始想用 log-distance 模型实时算但后来发现对这个项目的主要研究点来说标称 PL 和距离的耦合不是重点。先把路径损耗固定成一个值集中分析衰落统计和 IRS 增益等核心模型验证完再引入距离模型后面调参数会更轻松。3.2 信道样本生成与蒙特卡洛循环蒙特卡洛循环的核心是生成 K×M 个小尺度功率增益 W 和 K 个阴影功率因子 S。这里要注意阴影衰落有空间相关性同一架无人机到用户方向上不同反射单元的路径非常接近所以它们共享同一个阴影因子。而小尺度衰落各单元之间是独立变化的必须独立生成。这个“一条反射支路共享阴影、各单元独立小尺度”的处理方式是仿真中必须保持的物理一致性。rho sqrt(10^((Ptx_dBm - N0_dBm - PL_dB) / 10)); % 幅度缩放因子 gamma_th 10^(gamma_th_dB / 10); % 门限转换线性值 outage_cnt 0; snr_store zeros(1, Nmc); for trial 1:Nmc % 生成小尺度功率增益每个IRS单元一个Nakagami-m功率样本 W gamrnd(m, Omega/m, K, M); % 生成逆伽马阴影每架无人机共享一个阴影样本 S 1 ./ gamrnd(alpha, 1/beta, K, 1); % 计算等效信噪比 h_amp sqrt(S .* W) * rho; % S、W逐元素相乘再开方 total_amp sum(h_amp(:)); % 理想相位对齐后相干叠加 snr total_amp^2; snr_store(trial) snr; if snr gamma_th outage_cnt outage_cnt 1; end end p_out outage_cnt / Nmc;这个循环里有个很关键的小细节S .* W的维度和W不一样一个是 K×M一个是 K×1。MATLAB 会自动做扩展让每个阴影样本 S(k) 乘到这一架无人机所有 M 个反射单元的 W 上刚好实现了“同机同阴影异机异阴影”的物理设定。如果你的 MATLAB 版本较老或者你自己手写循环千万别把维度写错否则生成的信道就是一塌糊涂。3.3 智能反射面相移优化与多无人机合并智能反射面的“智能”二字在于每个反射单元可以调整自己的相位。数学上如果所有反射单元的信号在用户终端是同相叠加那么等效信道幅度就是幅度的代数和。我在代码里用sum(h_amp(:))直接完成这个操作这其实隐含了一个强假设所有反射链路的相位都被完美对齐了。实际实现中每架无人机 IRS 的相移是根据该无人机到用户的信道相位来设置的。真实信道会有相位噪声但仿真阶段先做理想相位对齐这样能给出性能上界。如果后面要研究非理想情况可以在每个反射链路上额外叠加一个随机相位误差 把相干求和改成sum(h_amp .* exp(1j * phase_err))再对合成信噪比取模。这是我后续准备扩展的方向之一。多无人机合并还有一个容易被忽略的点无人机群协同不等于把每架无人机的能量简单相加而是要关注不同无人机反射链路之间的阴影相关性。如果所有无人机飞得很近共享同一个大尺度阴影那么 S 就只生成一个而不是 K 个此时无人机多样性增益会明显下降。这正好对应了我后面要做的对比实验无人机群充分散开时的性能和聚在一起的性能差多少。4. 关键结果与性能分析4.1 中断概率随 m 值的变化趋势先看小尺度衰落参数 m 的影响。固定 alpha4、beta3、K4、M16扫描 m 从 0.8、1、2、3 到 5。仿真结果曲线整体上呈现一个合乎直觉的趋势m 越大中断概率单调下降。原因是 m 增大后每个反射链路的小尺度功率增益方差变小出现“深的衰落谷”的概率变低多个 IRS 单元相干叠加后的合成 SNR 更集中在高值区域。这里要特别说明的是m 从 0.8 增加到 1中断概率下降非常明显因为这是从“比瑞利更恶劣”到标准瑞利的跃迁。而 m 从 3 增加到 5曲线的改善幅度就放缓了。原因是大 m 值下信道已经接近静态此时的瓶颈不再是多径衰落而是阴影和路径损耗。这个饱和效应在无人机群场景下尤其明显因为它有 64 条反射路径4 架 × 16 单元做分集单独增加 m 对总体的边际效益会被阵列增益吞掉。实际调整参数时如果发现你的系统性能在某个 m 之上就不再改善就要想想是否值得花资源去改善信道质量。比如为无人机选择一条更稳定的飞行航线来提升 m如果提升有限倒不如增加无人机数量或者优化 IRS 单元数量投资回报比更高。4.2 阴影参数 α、β 对系统的影响逆伽马阴影的参数对系统性能的影响比 Nakagami-m 更“野”因为重尾特性会导致极端样本。我固定了 m2分别测了 alpha3 和 alpha6 两组同时让 beta 值跟着调整以保持阴影中值基本不变。结果显示alpha 越小中断概率曲线的“零散程度”越大也就是说单次仿真结果的抖动性很强。为什么会这样逆伽马分布尾巴衰减慢小 alpha 意味着有更大概率生成巨大的阴影功率 S这些异常样本对应“瞬时信道极差”的场景。在蒙特卡洛循环里这些极端样本虽然占比不大但会把平均 SNR 拉高同时把中断概率统计变得不稳定。这时候如果只跑 1e4 次蒙特卡洛结果很难收敛我最后把次数加到 1e5 甚至 5e5曲线才稳定下来。beta 参数的作用更偏向于缩放整体阴影水平。在其他条件不变时beta 增大相当于逆伽马分布向小值方向压缩也就是阴影效应减弱中断概率下降。有意思的是 beta 与 alpha 的耦合关系比较复杂因为它们共同决定分布的均值和方差。严格校准逆伽马阴影参数时最好先计算出分布的理论均值再反推 beta 值来匹配实际场景的中值功率而不是凭感觉乱给参数。4.3 无人机数量与 IRS 阵元数的组合增益这个锅最好用数据说话。我仿真了三种配置单无人机 16 单元、4 无人机各 16 单元、4 无人机各 64 单元。固定发射功率和门限统计中断概率。结果符合一个基本规律从单机到 4 机等效反射路径数变成原来的 4 倍再加上相干叠加信噪比增益远不止线性倍数。理论上完美相干叠加下总幅度与路径数成正比SNR 与路径数的平方成正比。从 16 条路径变成 64 条路径SNR 理论上提升 12 dB。当然这是在所有反射链路相位都对齐、路径损耗都相同的前提下得到的上界。但仿真结果也提示了一个陷阱K 增大到一定程度后中断概率下降速度变缓。原因有两方面一是阴影衰落引入了相关性如果多架无人机之间距离不够远它们的阴影样本并不完全独立分集增益会被抵消二是随着反射路径增多等效 SNR 上升到很高水平中断概率已经接近零继续增加无人机或阵元自然没有更多下降空间。所以纯堆数量不是最优解需要考虑覆盖范围和无人机间距。从工程上讲这给了一个很实际的设计思路如果你想用无人机群 IRS 覆盖一片区域与其每架无人机都塞几百个反射单元贵且重不如布置 6~8 架小型无人机携带 32~64 单元的小型 IRS通过合理分布获得分集和阵列增益。这个结论我是在仿真里先把 M 固定在 16、然后让 K 从 1 扫到 10 得到的K4 到 K6 的提升最大再往上就边际递减了。5. 实操经验与避坑指南5.1 MATLAB 仿真中的经典坑先说我踩过最深的坑随机数种子不固定。蒙特卡洛仿真结果波动大如果没有rng(2025)这样的种子锁定每次运行代码结果都不一样你就很难精确对比不同参数下的曲线。尤其是逆伽马阴影生成时极端样本偶尔出现可能直接改变整体结果。后来我习惯在每个仿真脚本开头先rng(shuffle)作为“正式跑”的设置但做调试和参数扫描时统一用固定种子这样才能保证可重复性。第二个坑是gamrnd的参数含义。MATLAB 的gamrnd(a,b)中 b 是尺度参数期望是 ab方差是 ab²。但我项目里用的逆伽马是从“速率参数”角度定义的所以必须先转换成1/beta作为尺度。如果直接把 beta 塞进去生成的阴影会偏得离谱。早期版本我跑出来的中断概率几乎接近 1排查了半天才发现是这里尺度参数搞反了。第三个坑是数组维度不匹配。我在 3.2 节里特意强调了S和W的维度关系一旦写成S 1./gamrnd(alpha, 1/beta, 1, K)行向量再和 K×M 的 W 做逐元素乘法时 MATLAB 会按广播规则生成 K×K 的矩阵结果完全错误。遇到这种问题我建议你把中间量打到工作区里用size()检查一下别裸奔。第四个坑是循环性能。蒙特卡洛次数一上去纯for循环跑起来非常慢。我的做法是尽量向量化在循环内只生成信道并计算 SNR如果是对某个参数扫描可以先用parforParallel Computing Toolbox并行跑不同参数点。但注意parfor里要保留rng的正确性最好在每个 worker 内用独立的随机流避免重复样本。5.2 正确性验证三件套做仿真最怕的是代码跑完不知道自己结果对不对。我总结了三个验证手段每次改模型都必须过一遍。第一是退化验证。把逆伽马阴影的alpha设得很大比如 1000同时把beta设成 alpha-1让阴影的均值收敛到 1、方差趋近 0此时阴影项退化为常数 1再把 Nakagami-m 的 m 设为 1此时小尺度衰落退化为瑞利。整个系统就退化为“纯瑞利多径 无阴影”的简单场景你可以把中断概率和理论瑞利单链路中断概率做对比。如果两者接近说明核心采样逻辑没大问题。第二是理论闭合解对比。在单反射单元、无 IRS 阵元增益的极限下复合信道的功率增益分布等于 Gamma×逆伽马我可以在论文里找类似模型的闭合表达式把理论中断概率和蒙特卡洛仿真曲线画在同一张图上重合度越高代码可信度越高。第三是收敛性测试。固定参数分别跑 Nmc 1e3、1e4、1e5、5e5观察中断概率的变化范围。如果从 1e4 到 1e5 的变化小于 5%基本可以认为结果收敛了。逆伽马阴影重尾严重收敛慢所以我在正式结果中至少跑 5e5 次并且每组参数用同一个随机种子保证相对比较的公平性。5.3 这个模型还能怎么扩展这个项目的底子搭好之后后面扩展空间非常大。我个人建议优先做“非理想相位”和“动态无人机位置”这两件事。非理想相位好理解把 IRS 相移优化从理想变成有误差可以模拟量化误差、时钟不同步、信道估计不准。具体做法是在每个反射单元的幅度上乘上exp(1j*phi_err)其中phi_err可以是均匀分布或高斯分布。这个改动只需要几行代码但能让结果更贴近实测。动态无人机位置更有意思。我现在是固定拓扑如果让每架无人机的位置在每次蒙特卡洛循环里按某种轨迹移动路径损耗和阴影参数也会跟着变。这个可以进一步研究“无人机协同波束对用户位置的敏感性”比如无人机间距拉大、高度变化对覆盖率的影响。甚至可以引入一个简单的优化循环搜索使中断概率最小的无人机三维坐标组合。不过这种优化通常比较费算力建议把信道生成函数改成可调用函数再用fminsearch或粒子群算法套一层。最后提一句如果你想把这份代码用于更进阶的研究建议把信道生成部分抽成独立函数输入参数 m、Omega、alpha、beta、K、M输出等效信道幅度。这样后面无论做功率分配、轨迹优化还是机器学习预测都能直接调用不用每次复制粘贴一大段脚本。我后来把代码重新整理成这种模块结构扩展起功能来确实省了很多事。