做认知无线电仿真这几年我有个越来越明显的体会能量检测器ED虽然在工程里用得最多但它并不是最靠得住的检测方案尤其在噪声不确定性明显的场景下性能衰减比想象中大得多。后来我把Pietra-Ricci指数检测器接进集中式数据融合协作频谱感知框架里在Matlab里跑通了一整套蒙特卡洛仿真效果比我原先用的能量检测方案稳定不少代码实现也意外地简洁。这篇文章就把我的实现思路、核心代码、门限标定方法和踩坑记录完整梳理一遍希望能给正在做协作频谱感知仿真的朋友提供一条可以直接参考的路径。1. 从单节点能量检测的软肋到协作感知的出现1.1 噪声不确定性与SNR Wall先聊聊我在实际仿真和文献阅读里反复遇到的一个问题能量检测器之所以简单是因为它只对接收信号做平方求和然后用这个能量值和噪声功率的门限比较。这里的关键在于门限必须精确知道当前环境下的噪声功率。可现实世界里的噪声从来不是安稳的——接收机热噪声本身会随温度漂移环境电磁干扰会带来短时波动前端滤波器带宽的变化也会改变噪声水平。一旦把噪声功率估计值算偏1dB能量检测器的虚警概率就可能从0.05飙到0.2或者反过来漏检概率居高不下。学术界把这种现象叫做SNR Wall意思是当噪声不确定度固定时无论你把采样时间拉多长都存在一个不可逾越的信噪比下限低于这个下限能量检测器根本无法可靠工作。这就像在嘈杂的餐厅里你明知道朋友会说某一句话但不知道他的音量用了哪种刻度来衡量结果再怎么集中精力也听不清。这个问题直接推动了两类研究方向一类是设计不依赖噪声功率的盲检测器比如基于协方差矩阵特征值的最大最小特征值检测器MME、协方差绝对值检测器CAV另一类是从单节点走向多节点协作通过多个认知用户的空间分集来对抗信道衰落和噪声干扰。Pietra-Ricci指数检测器其实正好站在这两类方向的交叉点上它用特征值构造统计量天然对噪声功率不敏感又能很方便地嵌入协作融合架构中。1.2 协作感知为什么能把这件事做成单节点频谱感知最大的瓶颈在于隐藏终端问题主用户信号可能正在发射但某个认知用户恰好处于深衰落区域它收不到信号就会误判频段空闲。多节点协作解决这个问题的思路很直接多个用户分布在不同的地理位置同一个主用户信号经历几条独立衰落的信道只要不是所有信道同时深衰落总有人能可靠地感知到信号。协作架构大体分两类。一类是分布式协作节点之间互换本地判决或统计量通过共识算法达成全局判决另一类是集中式协作所有节点把数据上报给一个融合中心由融合中心完成最终判决。本文选的集中式方案从Matlab仿真角度来说更直观从理论分析上也更好下手。融合中心相当于把所有节点看成了一个虚拟的大孔径接收机数据融合做得好分集增益非常明显。2. Pietra-Ricci指数检测器的原理推导2.1 从协方差矩阵到特征值分布Pietra-Ricci指数原本是统计学里衡量两个概率分布之间差异的度量后来被引入到信号检测领域用来刻画协方差矩阵的结构变化。先交代一下为什么协方差矩阵能用来做频谱感知。假设系统里有K个协作节点每个节点单天线融合中心在一个感知周期内拿到K个节点、M个快拍的接收数据排列成K×M的矩阵X。样本协方差矩阵就是[ \mathbf{R}x \frac{1}{M}\sum{m1}^{M}\mathbf{x}(m)\mathbf{x}^H(m) \frac{1}{M}\mathbf{X}\mathbf{X}^H ]在主用户不存在的H0假设下(\mathbf{x}(m))只包含噪声如果噪声是空时白噪声那么(\mathbf{R}_x)在统计意义上趋近于(\sigma_n^2\mathbf{I})是个对角阵所有特征值都落在(\sigma_n^2)附近。一旦主用户信号出现H1假设各节点收到的信号来自同一个主用户源节点间样本存在相关性(\mathbf{R}_x)不再是对角阵其特征值结构也会发生明显变化少数几个大特征值对应信号子空间其余特征值对应噪声子空间整体上特征值会展布开而不是挤在一起。所以检测主用户是否存在这个问题可以在数学上转化为判断协方差矩阵是否偏离白噪声矩阵结构。而特征值分布的发散程度就是这个偏离程度的天然度量。Pietra-Ricci指数检测器做的就是把这种特征值发散程度量化成一个标量统计量。2.2 PR检测统计量的构造与判决准则我采用的PR检测统计量定义如下设样本协方差矩阵(\mathbf{R}x)的K个特征值为(\lambda_1, \lambda_2, ..., \lambda_K)令(\bar{\lambda} \frac{1}{K}\sum{i1}^{K}\lambda_i)则PR指数为[ T_{\rm PR} \frac{2}{K} \sum_{i1}^{K} \frac{(\lambda_i - \bar{\lambda})^2}{(\lambda_i \bar{\lambda})^2} ]这个式子的物理含义很清晰。如果所有特征值都相等也就是(\lambda_i \bar{\lambda})对每个i都成立那么分子为0(T_{\rm PR}0)对应纯噪声场景。如果特征值差异很大比值项就会变大(T_{\rm PR})随之增大说明协方差矩阵明显偏离白噪声结构支持主用户存在的判决。要注意这个统计量有两个宝贵的性质。第一它是尺度不变的把所有特征值同时放大或缩小一个常数倍分子分母同阶缩放比值不变。这意味着门限不需要根据噪声功率调整理论上只要噪声还是白噪声哪怕绝对功率变了统计量的分布都不变。第二它不涉及矩阵求逆只做特征值运算数值上比某些需要估计噪声功率或求逆的检测器稳定得多。判决准则就是经典的二元假设检验预先通过蒙特卡洛在H0下标定门限(\gamma)当(T_{\rm PR} \gamma)判为H1主用户存在否则判为H0主用户不存在。门限的选择对应虚警概率(P_f P(T_{\rm PR} \gamma | H_0))。2.3 与能量检测、MME检测器的对比为了说明为什么我会在协作感知场景里选中PR检测器这里列出三种检测器的实测感受对比检测器统计量依赖对噪声功率的敏感度计算复杂度我的评价能量检测ED信号能量与噪声功率之比高噪声不确定时性能崩塌最低O(KM)大快好省但在低信噪比和噪声不确定场景不敢用最大最小特征值比MME(\lambda_{\max}/\lambda_{\min})不敏感无需噪声功率中需特征值分解经典方案但对个别最大特征值波动敏感Pietra-Ricci指数PR所有特征值与均值偏离程度不敏感尺度不变中需特征值分解全局刻画特征值展布抗单特征值突变融合场景更稳MME只用了最大和最小两个特征值意味着它主要关注特征值分布的两个极端。在实际小样本条件下最大特征值的波动往往很大虚警概率会偏高。PR指数把所有特征值都纳入统计量相当于对特征值分布做了全局衡量遇到极端特征值跳动时不会那么敏感。协作节点数增多以后PR指数的特征值维度也增高这种全局统计量的优势会更明显。3. 集中式数据融合架构与仿真系统模型3.1 系统模型K个节点、M个快拍、两种假设先明确仿真里用的系统模型。设有K个协作认知节点每个节点配置单根天线所有节点将采样数据上报到一个融合中心。融合中心在一个感知时隙内采集M个快拍。在H0假设下第k个节点的第m个采样为[ x_k(m) n_k(m) ]在H1假设下[ x_k(m) h_k s(m) n_k(m) ]其中(s(m))是主用户信号建模为零均值复高斯随机变量功率为(\sigma_s^2)(n_k(m))是零均值复高斯噪声功率为(\sigma_n^2)(h_k)是主用户到第k个节点的信道增益。如果按瑞利衰落建模(h_k)是循环对称复高斯随机变量即实部和虚部独立同分布、各占一半功率如果按确定性信道建模可以让(h_k1)验证算法基线性能。信噪比定义为接收端单个节点上的平均信噪比[ \mathrm{SNR} \frac{\mathbb{E}[|h_k|^2]\sigma_s^2}{\sigma_n^2} ]每个节点单天线、K个节点协作融合中心拿到的是K×M的复基带数据矩阵(\mathbf{X})。假设主用户信号带宽远小于中心频率可以采用窄带平坦衰落模型所有节点的信号包络在M个快拍内保持不变或独立变化都行仿真中按块衰落处理。3.2 融合中心的软融合策略与全局判决流程集中式数据融合可以分成硬融合和软融合两个层次。硬融合是每个节点先做本地判决只把0/1结果发给融合中心融合中心再按或、与或K-out-of-N规则合并。这种方案传输开销小但在本地判决阶段就已经损失了信息尤其是低信噪比环境下的细微相关性单节点根本判不出来。软融合则是节点把采样数据或足够的统计量上传给融合中心由中心做全局最优判决。我用的是软融合中的协方差域融合节点上传原始采样序列融合中心直接构造全局样本协方差矩阵。之所以在仿真里采用原始数据上传是因为它最直观也最能体现数据融合四个字的含义。实际工程里可以退化为节点本地算部分相关矩阵再上传但那是另外一套权衡了。融合中心拿到矩阵(\mathbf{X})后按下面流程执行根据K×M矩阵(\mathbf{X})计算K×K样本协方差矩阵(\mathbf{R}_x \mathbf{X}\mathbf{X}^H / M)对(\mathbf{R}_x)做特征值分解得到特征值(\lambda_1, ..., \lambda_K)代入PR指数公式计算(T_{\rm PR})将(T_{\rm PR})与门限(\gamma)比较输出全局判决结果。这个流程把K个节点在空间上综合成了一个虚拟天线阵列节点间的相关性被协方差矩阵完整保留。这正是集中式融合相对单节点检测的核心优势。4. Matlab实现模块、代码与参数配置4.1 顶层仿真结构与参数初始化整个Matlab仿真我拆成了三个模块参数初始化与门限标定、信号生成与协方差计算、性能统计与绘图。先给出一份可以直接运行的顶层脚本框架所有随机实验都建议开头设置随机种子保证可复现。%% 集中式协作频谱感知 - Pietra-Ricci指数检测器 % 仿真设置 clear; clc; close all; rng(42); % 随机种子确保可复现 K 8; % 协作节点数 M 2000; % 每个节点采样快拍数 Pf 0.05; % 目标虚警概率 num_mc 5000; % 蒙特卡洛次数 snr_dB -20:2:0; % 信噪比扫描范围 % 生成H0下的检测统计量用于标定门限 T_h0 zeros(num_mc, 1); for mc 1:num_mc X sqrt(0.5) * (randn(K, M) 1i*randn(K, M)); % 纯复高斯噪声 R (X * X) / M; lam eig(R); T_h0(mc) PR_index(lam); end gamma quantile(T_h0, 1 - Pf); % 经验门限 fprintf(Pf %.3f 对应的门限 gamma %.4f\n, Pf, gamma);参数里K8的意思是8个协作节点M2000表示每个节点上报2000个采样点。在认知无线电场景里频谱感知时隙通常是毫秒到几十毫秒量级采样率几MHz的话M取2000到5000是比较常见的范围。K越大空间分集越强但节点间数据同步和上报开销也越大。4.2 门限标定与检测概率计算门限标定是检测器仿真最容易出错的地方。PR指数是尺度不变的门限可以在噪声功率归一为1的条件下标定之后无论仿真中噪声功率怎么变门限都直接复用。虚警概率的统计意义是在纯噪声条件下统计量超过门限的概率。所以我用5000次蒙特卡洛生成H0分布取它的(1-Pf)分位数作为门限。%% 计算检测概率 Pd zeros(length(snr_dB), 1); for i 1:length(snr_dB) snr db2pow(snr_dB(i)); sig sqrt(snr * 0.5); nois sqrt(0.5); count 0; for mc 1:num_mc % 瑞利衰落信道 h (randn(K, 1) 1i*randn(K, 1)) / sqrt(2); % 主用户信号零均值复高斯功率 snr s sig * (randn(1, M) 1i*randn(1, M)); % 噪声功率 1 N nois * (randn(K, M) 1i*randn(K, M)); % 接收数据 X h * s N; R (X * X) / M; lam eig(R); T_PR PR_index(lam); if T_PR gamma count count 1; end end Pd(i) count / num_mc; end代码里signal的生成直接用sqrt(snr*0.5)*(randn1i*randn)这保证信号样本的功率期望是snr。噪声类似地保证功率期望为1。这样SNR的定义就是线性的功率比不需要额外归一化。特别提一句很多时候新手直接在Matlab里用randn生成所谓噪声不乘sqrt(0.5)最终会导致信号噪声功率比偏离预期曲线整体平移几个dB。4.3 PR指数计算函数与核心细节PR指数计算函数只有几行但里面有一个数值细节值得单独说明。function T PR_index(lambda) % lambda: 样本协方差矩阵的特征值向量 % 返回: Pietra-Ricci指数检测统计量 lambda lambda(:); lambda max(lambda, eps); % 防止极小特征值导致数值问题 N length(lambda); mu mean(lambda); T 2 * sum((lambda - mu).^2 ./ (lambda mu).^2) / N; endmax(lambda, eps)这行是我在实际调试里加上的。样本协方差矩阵理论上半正定特征值非负但浮点运算可能算出接近0甚至微小负数的特征值。如果不做下限保护分母里出现负数或零PR指数可能出现无穷大或NaN。对于K8的小矩阵这种概率不高但蒙特卡洛跑上万次之后总会撞上所以建议保留这个保护。特征值分解我用的是eig(R)。对于Hermitian矩阵也可以改用svd(R)取奇异值的平方不过eig专门针对对称矩阵做了优化速度更快内存占用也更低。如果K超过50可以考虑用eig的vector选项只求特征值进一步减少计算量。4.4 代码中容易被忽略的细节我刚跑通这个仿真时踩了几个坑这里集中提醒第一信号和噪声都要用复基带模型。认知无线电接收机做数字信号处理时下变频后是I/Q两路所以仿真里的数据必须是复信号。如果用实数高斯模拟协方差矩阵的秩特性和统计分布会和实际接收机对不上检测性能自然也会失真。第二瑞利衰落信道的功率归一化。h (randn(K,1)1i*randn(K,1))/sqrt(2)这个写法不是可有可无/sqrt(2)保证了(|h_k|^2)的期望为1。如果漏掉这一步实际信噪比会偏小3dB最终得到的检测概率曲线会整体右移你还会误以为是算法性能差。第三蒙特卡洛次数不是越大越好但5000次是底线。尤其是在标定门限时Pf0.05意味着H0分布尾部只有5%的概率5000次里只有约250个样本在尾部门限估计的波动已经比较明显。想要更平滑的曲线建议门限标定用10000次检测概率统计用5000到10000次。我的经验是如果门限标定次数不够低信噪比段的虚警概率会明显偏离0.05看起来就像检测器失效了。5. 仿真结果检测概率、节点数与噪声稳健性5.1 不同SNR下的检测性能曲线用上面参数跑出来的结果趋势是非常典型的。当SNR从-20dB逐步上升到0dB时检测概率Pd会从接近0的位置开始爬升最终到达1。在K8、M2000、Pf0.05的参数下大约在-12dB附近Pd就能达到0.9以上。这里要解释一下为什么低信噪比段PR检测器表现不错。8个协作节点让协方差矩阵的维数达到8信号子空间虽然只有一个主特征值但它的出现会把整个特征值分布的失衡程度拉大。PR指数把所有特征值相对均值的偏差都统计进来因此微弱信号导致的微小失衡也能被累积放大。相比之下单节点能量检测在-12dB时需要非常精确的噪声功率估计实测里很难做到。当然如果M很小比如只有200个快拍样本协方差矩阵的估计误差会变大PR指数的分布会展宽门限标定和检测性能都会下降。这是所有基于协方差矩阵的检测器的共同特点小样本条件下特征值波动大统计量分布重叠严重。5.2 协作节点数K如何影响性能我对比过K2、4、8、12这四种配置固定M2000、Pf0.05。结论是协作节点数从2增加到8时检测概率曲线的增益非常明显大约有4~6dB的SNR增益但从8增加到12时增益开始变缓。这个现象背后的物理原因很直观空间分集增益随节点数增加而增加但边际增益是递减的。另一个容易被忽略的点是K越大协方差矩阵维数越高同样的M下协方差估计质量会下降。换句话说如果节点数翻倍但每个节点的采样快拍数不变特征值散布会变大H0分布尾部更肥门限也会抬高反而削弱检测能力。所以实际工程中K和M需要联动权衡并不是节点越多越好。我从仿真里得到的经验是K8、M2000左右是一个性能和开销比较平衡的工作点。5.3 抗噪声不确定性验证这一节是我个人最看重的验证。做法很简单在H1检测时不再假定融合中心知道准确的噪声功率而是给噪声样本乘一个随机的功率偏差因子。% 噪声不确定度验证每个蒙特卡洛循环随机化噪声功率 noise_power_factor db2pow(randn * 1); % 标准差1dB的随机偏差 N sqrt(0.5 * noise_power_factor) * (randn(K, M) 1i*randn(K, M));这里的randn*1表示噪声功率在0dB附近做±1dB左右的随机扰动。对能量检测器来说这种扰动足以让虚警概率和检测概率同时失控但对PR检测器来说因为统计量本身是尺度不变的噪声功率的绝对变化不会改变特征值之间的比值结构检测性能几乎不受影响。我在相同条件下对比过PR检测器加不加噪声不确定度Pd曲线差距小于1个百分点这在实际检测问题里是个非常宝贵的性质。6. 工程实现中的几个坑与性能优化建议6.1 协方差矩阵的病态与正则化处理当节点数K接近甚至超过快拍数M时样本协方差矩阵会变成秩亏矩阵出现大量零特征值。这种情况下PR指数的分母会出现问题统计量会整体偏向很大的值虚警概率飙升。解决思路是正则化在计算协方差矩阵时加一个很小的对角加载项也就是用(\mathbf{R}_x \epsilon\mathbf{I})代替(\mathbf{R}_x)。(\epsilon)的选择有个经验法则取(\epsilon 10^{-5}\cdot\mathrm{trace}(\mathbf{R}_x)/K)既能保持特征值结构不变又能避免秩亏问题。如果M远大于K正则化就不是必须的。所以仿真前先看一下K×M的维度关系别急着加对角加载也别完全不处理。6.2 蒙特卡洛次数与门限精度的取舍门限标定的本质是用经验分位数估计理论分位数样本量不足时误差很大。我在实际测试中发现一个现象如果只做2000次蒙特卡洛标定门限运行完整仿真后测得的实际虚警概率可能是0.08而不是目标的0.05让人误以为检测器有问题。增大到8000到10000次后实际虚警概率才稳定收敛到0.05附近。另外低信噪比段检测概率的蒙特卡洛方差也更大因为此时检测统计量恰好落在门限附近的概率很高每次仿真掷硬币的结果对Pd影响很大。我的建议是先跑一个小规模实验确定大概曲线形状再在关键SNR点比如Pd在0.5附近的点加大蒙特卡洛次数得到更精确的估计。6.3 从仿真到实机部署的扩展方向最后聊一下代码往实机方向走时需要考虑的扩展。首先是数据上报压缩原始采样上传带宽很大实机里通常改为每个节点本地计算采样自相关并上传但要注意只上传自相关会丢失节点间互相关信息。一个折中方案是每个节点同时上传自身数据与少数关键邻居节点的互相关估计这本质上是分布式协方差估计问题。其次是信道模型的扩展窄带平坦衰落模型适合窄带频谱感知如果是宽带信号需要把数据分成多个子带分别做PR指数检测再在各子带间融合。这个扩展在Matlab里并不难只要把顶层循环改成频点循环即可。还有一个我最近在尝试的方向是把PR指数和能量检测器做成自适应选择先估计噪声环境的不确定度不确定度低就用能量检测计算量小不确定度高就切到PR检测器鲁棒性强。两个检测器共享同一份协方差矩阵特征值PR指数代码可以直接复用几乎不增加额外开销。这种混合方案我觉得在实际系统里会比单用任何一种检测器都更实用。