基于Pietra-Ricci指数的协作频谱感知Matlab仿真：集中式融合实现

协作频谱感知这个方向我在认知无线电的仿真项目里反复接触过最头疼的永远不是怎么判断主用户在不在而是噪声功率稍微变一变传统能量检测就彻底失灵。这次分享的项目是一个挺有代表性的方案用Pietra-Ricci指数检测器做集中式数据融合的协作频谱感知全部用Matlab实现。它解决的核心问题是——在没有噪声功率先验、没有主用户信号波形先验的情况下如何靠多个节点的数据融合把弱信号检测出来。这个项目比较适合正在做频谱感知仿真、认知无线电课题或者毕业设计选相关方向的同学。Pietra-Ricci简称PR指数本质上是一个衡量两个概率分布差异有多大的指标把它搬进频谱感知里就是通过比较只有噪声时特征值分布长什么样和信号加噪声时特征值分布长什么样来判断主用户是否存在。比起能量检测和匹配滤波这类盲检测器不依赖噪声功率精确已知实际信道环境下更抗造。下面我从设计思路、数学原理、Matlab实现到常见坑完整拆一遍。1. 项目整体设计与思路拆解1.1 协作频谱感知为什么需要集中式数据融合先说说为什么单节点频谱感知不够用。单个感知节点在真实环境里会遇到三个致命问题阴影衰落把主用户信号压到噪声底以下、多径衰落造成深度谷点、还有隐藏终端问题——节点位置正好在覆盖盲区主用户信号根本到不了它那里。这时候节点自己判断信道空闲实际上是在跟主用户抢频段干扰就产生了。协作感知的思路很简单多放几个节点在空间不同位置各自感知后再把信息汇总利用空间分集来对冲个别节点的糟糕信道条件。只要不是所有节点同时处于深度衰落整体判决就有很大概率做对。集中式数据融合是这里面最经典的一种架构所有次用户SU把本地感知结果上传到一个融合中心FC由FC做最终判决。它跟分布式协作的区别在于FC拥有全局信息理论上能达到最优检测性能算法设计、性能分析都更可控。代价是FC单点故障风险和回传链路开销但这个代价在理论研究和性能基准验证阶段完全可以接受。我这次实现的就是这种多节点采样-汇聚-中心判决的链路。1.2 为什么选Pietra-Ricci指数检测器而不是能量检测能量检测为什么在工程里不好用因为它需要知道噪声功率才能设定正确的判决门限。接收信号能量 信号能量 噪声能量只有把门限设在噪声能量之上、信号能量之下才能分得开。问题是噪声功率不是恒定的——温度变化、射频前端增益漂移、邻频干扰都会让噪声功率在±1dB甚至更大范围内波动。门限设低了虚警爆炸设高了漏检严重。匹配滤波检测性能好但它需要知道主用户信号的完整先验波形、调制方式、定时这在非合作的频谱感知场景里基本不现实。循环平稳检测不需要先验但要积累大量样本计算复杂度高实时性差。PR指数检测器走的是另一条路把信号检测问题转成分布一致性检验问题。它不需要知道噪声功率的绝对值也不需要信号波形只需要一个理论参考——纯噪声环境下接收信号协方差矩阵特征值的分布。这个参考分布是可以用数学推导或者离线仿真精确摸清的。检测时把实测特征值分布和参考分布一比较差异大就判H1差异小就判H0。噪声功率的绝对大小不会改变特征值分布的形状归一化之后所以天然免疫噪声不确定性。1.3 集中式数据融合下的PR检测整体处理链路这套系统的完整信号流是这样的每个感知节点在感知时隙内采集一段复基带信号本地做简单的预处理去直流、归一化然后把原始采样数据或者本地采样协方差矩阵上传给融合中心。融合中心拿到所有节点的数据后拼接成一个全局数据矩阵计算全局采样协方差矩阵做特征值分解再算特征值经验分布和理论噪声分布的PR距离。最后把这个距离跟预设门限比较超过门限判主用户存在否则判信道空闲。2. Pietra-Ricci指数的数学原理与检测门限设计2.1 PR指数的定义和直观理解Pietra-Ricci指数是个度量两个概率分布之间差异的指标。给定两个累积分布函数F(x)和G(x)PR距离定义为[ D_{PR}(F, G) \frac{1}{2} \int_{-\infty}^{\infty} |F(x) - G(x)| dx ]乘1/2是为了让距离范围落在[0, 1]区间。两个分布完全重合时FG距离为0完全分开时距离趋近1。它衡量的是两条CDF曲线之间夹着的面积。这里要区分一下PR距离和KS统计量。KS统计量取的是两条CDF的最大垂直距离sup|F-G|只关注哪里差得最狠PR距离看的是整体差了多少。在频谱感知场景里噪声不确定性带来的特征值分布变化往往是整体性的展宽和平移而不是某一单点的陡变所以PR距离对这类变化的敏感度更好统计量更平滑、更抗单点异常特征值的干扰。2.2 检测统计量怎么从数据里算出来先说接收信号模型。H0假设下第i个节点接收到的复基带信号是纯噪声[ x_i(n) w_i(n), \quad w_i(n) \sim \mathcal{CN}(0, \sigma^2) ]H1假设下是主用户信号加噪声[ x_i(n) s_i(n) w_i(n) ]每个节点采集L个快拍。集中式融合这里有两种做法一种是节点直接上报原始数据FC把所有节点的数据竖着拼起来得到一个(M·N)行、L列的矩阵X另一种是每个节点本地先算协方差矩阵FC把所有协方差矩阵平均。第一种更接近数据级融合信息损失最小我实现的就是这种。FC拿到全局数据矩阵后计算全局采样协方差矩阵[ R \frac{1}{L} X X^H ]对R做特征值分解得到一组特征值λ1, λ2, ..., λp。然后把特征值按升序排列再除以特征值均值做归一化[ \tilde{\lambda}i \frac{\lambda_i}{\frac{1}{p}\sum{j1}^p \lambda_j} ]这个归一化是盲检测的关键。噪声功率σ²对每个特征值的影响近似是等比例的归一化之后纯噪声下的特征值分布形状就跟σ²无关了只跟节点数、快拍数、协方差矩阵维度有关。接下来构造经验CDF。把归一化后的特征值(\tilde{\lambda}_i)作为横轴经验CDF的纵轴值取((i-0.5)/p)就得到了实测特征值分布(\hat{F}(x))。参考分布F0(x)怎么来理论上有两种途径。一是用Marchenko-PasturMP律它给出白噪声协方差矩阵特征值的渐近分布二是用离线仿真模板——在纯噪声假设下用同一套节点数、快拍数参数跑大量蒙特卡洛把特征值经验CDF存下来做模板。工程上我强烈推荐第二种因为MP律是渐近结果在小样本、维度不高时偏差不小而仿真模板天然适配你的实际参数谁用谁知道。最后检测统计量就是实测经验CDF和参考CDF之间的PR距离[ T \frac{1}{2} \int |\hat{F}(x) - F_0(x)| dx ]Matlab里用trapz做梯形积分就能算出来。2.3 检测门限与虚警概率控制PR统计量的分布没有闭合解析表达式门限必须靠离线蒙特卡洛标定。做法是在纯噪声假设下按你的场景参数节点数M、每节点快拍数L、协方差矩阵维度p生成大量H0数据每个样本算出一个PR统计量T。把这些T的分布累积起来取它的(1-Pfa)分位数就是你要的判决门限γ。[ \gamma F_T^{-1}(1 - P_{fa}) ]比如要Pfa0.1就把纯噪声下所有T值排序取第90百分位数。这个标定过程要注意两点一是蒙特卡洛次数至少20000次起步否则高分位数抖动很大二是标定用的随机数流和后面性能仿真用的随机数流要分开不然会引入乐观偏差。门限γ只跟系统参数M、L、p有关跟SNR无关所以可以预先离线算好线上检测时直接查表。实测下来在节点数M4、每节点快拍数L1024的配置下Pfa0.1对应的门限大概在0.3到0.5这个量级具体以你的离线标定结果为准。3. Matlab仿真系统实现全流程3.1 仿真场景与参数设置先把仿真参数列清楚。主用户信号用复正弦或者QPSK调制信号都行关键是信号要经过信道——我建议至少加一个频率平坦衰落信道不然太干净了体现不出检测器的优势。参考参数配置% 仿真参数配置 M 4; % 协作节点数 N 8; % 每个节点的接收天线数/协方差矩阵维度 L 1024; % 每个节点采集的快拍数 SNR_dB -15:2:5; % 信噪比扫描范围 Pfa_target 0.1; % 目标虚警概率 num_mc 5000; % 性能仿真的蒙特卡洛次数 num_th 20000; % 门限标定的蒙特卡洛次数 rng(2024); % 固定随机种子保证实验可复现这里要解释一下N的选取。N是指参与特征值分解的协方差矩阵维度在集中式融合里它等于所有节点数据堆叠后的总行数。如果每个节点是单天线那N就等于节点数M乘以每节点输出信号的路数。N太小特征值个数太少经验CDF分不出精细形状N太大协方差矩阵估计需要的快拍数L也得跟着涨否则矩阵不满秩。经验公式是L至少要是N的4到5倍我一般取L8N以上。3.2 节点本地处理与融合中心判决每个节点的任务是采集数据、简单预处理、上传。这里有个细节节点上传的是I/Q复数据不是能量值这样才能在FC做数据级融合。数据量确实大但这是性能上限的参考方案。FC的处理分四步第一步拼接数据矩阵。假设第i个节点上传的是一个N_sub×L的矩阵Y_iFC把所有节点的数据垂直拼接得到全局矩阵X尺寸是(M·N_sub)×L。第二步计算全局采样协方差矩阵并做特征分解R (X * X) / L; lambda eig(R); lambda sort(real(lambda), ascend);注意eig对复矩阵可能返回复数特征值实际上协方差矩阵是Hermitian半正定的特征值一定是实数取real是保险动作。第三步归一化特征值并计算经验CDF和PR统计量lambda lambda / mean(lambda); p length(lambda); Fhat ((1:p) - 0.5) / p; % 参考分布用离线模板这里直接调用 F0_ref reference_cdf(); % 长度为p的列向量预先离线生成 T 0.5 * trapz(lambda, abs(Fhat - F0_ref));第四步与门限比较if T gamma_threshold decision 1; % 判主用户存在 else decision 0; % 判信道空闲 end3.3 参考CDF模板的离线生成这一步是整个检测器能不能work的关键我单独拿出来讲。参考CDF模板必须在H0假设下生成也就是严格纯噪声不含任何信号分量。生成方法和线上流程一模一样只是把特征值换成纯噪声下的特征值function F0 generate_reference_cdf(M, N_sub, L, num_sim) p M * N_sub; F0_accum zeros(p, 1); for k 1:num_sim X (randn(p, L) 1j * randn(p, L)) / sqrt(2); R (X * X) / L; lambda eig(R); lambda sort(real(lambda), ascend); lambda lambda / mean(lambda); % 把每条H0样本的经验CDF值累加 F0_accum F0_accum ((1:p) - 0.5) / p; end F0 F0_accum / num_sim; % 同时把对应的横轴lambda也保存下来 end这里有个实操细节CDF模板的横轴是归一化特征值但每次蒙特卡洛生成的特征值都不完全一样不能直接对CDF值做平均要先在公共横轴上做插值再平均。更简单的做法是只在特征值位置上计算PR距离参考CDF用MP律的理论值或者像我实际项目里那样把纯噪声的特征值分布用核密度估计拟合成一条光滑曲线线上检测时用这条曲线插值。总之目标是得到一条稳定的纯噪声特征值CDF曲线。3.4 蒙特卡洛性能仿真主循环门限标定好、参考CDF模板准备好之后性能仿真就简单了Pd zeros(length(SNR_dB), 1); Pfa_sim zeros(1, 1); % 先跑H0验证虚警是否落在目标值附近 for mc 1:num_th X generate_noise_only(M, N_sub, L); T compute_pr_statistic(X); T_h0(mc) T; end gamma quantile(T_h0, 1 - Pfa_target); Pfa_sim mean(T_h0 gamma); fprintf(标定Pfa%.4f, 实测Pfa%.4f\n, Pfa_target, Pfa_sim); % 再跑H1得到不同SNR下的检测概率 for idx 1:length(SNR_dB) snr SNR_dB(idx); det_count 0; for mc 1:num_mc X generate_signal_plus_noise(M, N_sub, L, snr); T compute_pr_statistic(X); if T gamma det_count det_count 1; end end Pd(idx) det_count / num_mc; end画图部分用semilogy还是plot看个人偏好我习惯画两条曲线一条是Pd vs SNR固定Pfa0.1另一条是ROC曲线固定SNR扫描门限。ROC曲线扫描门限时要用同一批H0和H1的统计量只改判决门限不要重新生成数据。4. 典型实验结果与性能分析4.1 节点数对检测性能的影响我在固定每节点快拍数L1024、目标Pfa0.1的条件下分别跑了M1、2、4、8四个配置。结果符合预期单节点在SNR-10dB左右检测概率就开始明显下滑M4时这个拐点能往左移3到5dBM8比M4又提升1到2dB但边际收益明显递减。这个现象背后的道理是协作节点数增加协方差矩阵维度p变大特征值个数变多经验CDF的形状更稳定、噪声波动被平均掉所以PR距离的H0分布更集中同一门限下H1分布更容易分离。但超过一定数量后新增节点提供的空间分集增益趋于饱和反而因为需要的快拍数更多、回传开销更大投入产出比下降。实际工程里M4到6是个比较划算的区间这也是很多文献里爱用M4的原因。4.2 快拍数L对检测性能的影响把M固定在4L从256扫到4096结果就是检测概率随L增大全面提升尤其在低SNR区提升幅度非常明显。-15dB SNR下L256时基本测不到L4096时检测概率能到0.7以上。原因不复杂快拍数决定协方差矩阵估计的准确度。L越大采样协方差矩阵越接近真实协方差矩阵特征值分布越稳定经验CDF和参考CDF在H0下的偏差越小统计量方差越小检测器分辨率越高。这也是PR检测器的一个特性——它在低SNR下换取性能的方式就是堆样本而且堆样本的效果比能量检测更显著因为它利用的是特征值分布的集合效应不只是能量总量。4.3 与能量检测的对比噪声不确定性场景表格对比一下更直观场景能量检测噪声功率精确已知能量检测噪声功率有±1dB误差PR指数检测器检测概率SNR-12dB, Pfa0.1约0.85约0.40约0.78是否需要噪声功率先验需要需要但不可靠不需要是否需要信号波形先验不需要不需要不需要计算复杂度极低极低中等特征值分解这个表里的具体数值是某一组参数下的实测结果换个参数绝对值会变但趋势是稳定复现的噪声功率精确已知时能量检测略胜一筹一旦噪声功率估计有偏差能量检测性能暴跌PR检测器几乎不受影响。这就是PR检测器存在的最大价值——在非理想信道环境下做稳健检测。5. 常见问题与排查实录5.1 特征值出现NaN或Inf这个坑我踩过不止一次。根因基本是两个一是数据没归一化直接拿原始ADC采样值算协方差矩阵动态范围一大就容易溢出二是协方差矩阵不满秩——快拍数L小于协方差矩阵维度p导致R是奇异的特征值里出现0或负得很离谱的值。解决方法是先对每个节点的数据做均值去除和方差归一化再上传同时严格要求L ≥ 8p在代码里加个断言assert(L 8 * p, 快拍数过少协方差矩阵可能不满秩);如果还是出现NaN检查数据里有没有NaN源头比如randn生成时没设种子导致复现异常或者信道系数里不小心生成了Inf。5.2 实测虚警概率和目标Pfa对不上这是门限标定环节最常遇到的问题。我遇到过实测Pfa比目标Pfa翻倍的情况排查后发现是参考CDF模板和线上检测用的参数不一致——标定门限时用的节点数是4线上检测时实际节点数变成了6统计量分布整体平移了门限自然失效。另一个原因是蒙特卡洛次数太少导致门限不精确。纯噪声下PR统计量的分布拖尾比较长2000次蒙特卡洛去估0.01分位数误差能有30%以上。解决办法门限标定至少20000次分位数越极端次数要求越高。另外检查一下你在标定门限和实际检测时是不是都用了相同的特征值归一化步骤。特征值忘记除以均值统计量会整体变大好几倍门限直接作废。5.3 融合方式到底选硬融合还是软融合有些入门读者会问为什么不直接用多数表决这种硬融合——每个节点本地判决FC统计一下有几个判H1超过一半就判占用。硬融合的好处是回传开销极小只传1bit但代价是损失了软信息。一个节点SNR很高判得很有把握跟一个节点SNR刚好在临界点、碰巧判对在硬融合里权重一样这显然不是最优。本项目用的数据级融合上传原始样本FC构造全局协方差矩阵是性能上限但实现成本也最高。中间路线是软融合每个节点上传自己的PR统计量FC对统计量做加权合并。我给个建议排序做理论基准研究用数据级融合做工程仿真选软融合做低功耗实现才考虑硬融合。5.4 实操心得验证检测器正确性的快速方法最后分享一个我自己的调试套路。新写好的PR检测器先别急着上蒙特卡洛用三个快速测试验证正确性第一步纯噪声下跑500次统计T的均值应该是一个很小的值比如0.05以下如果发现均值比门限还大说明参考CDF模板或者归一化逻辑有bug。第二步高SNR比如10dB下跑100次T均值应该明显大于门限如果H1下的统计量反而比H0还小检查特征值排序方向、CDF计算有没有搞反。第三步用同一批数据分别跑PR检测器和能量检测在SNR中等0dB时两者的判决应该高度一致只有噪声不确定场景下PR才表现出优势。这三个测试通过了再跑完整的蒙特卡洛曲线基本不会出大问题。我个人在实际项目里体会最深的一点是PR检测器的优势不是靠更高端的数学赢来的而是靠换了一个更稳健的视角赢来的。能量检测盯着能量的绝对大小PR检测器盯着分布的形状差异后者对环境的未知因素天然免疫。这个思路不止适用于频谱感知任何信号检测问题——瞬态信号检测、异常检测、故障诊断——只要有相对干净的参考分布都可以想想能不能用上PR距离这个工具。扩展方向上加权软融合、跟深度特征结合做非高斯噪声下的检测都是可以继续挖的方向。做仿真时把门限标定和性能评估分成两套独立流程养成这个习惯你的结果会可靠很多。