做认知无线电仿真的朋友十有八九都跟频谱感知打过交道。如果把“Pietra-Ricci指数检测器”和“集中式数据融合”这两个词放一起你可能觉得很冷门——但它其实是一条非常实用的实现路径多个次用户把本地能量统计量汇聚到融合中心融合中心不需要复杂的协同算法只用一个分布距离指标就能完成判决整条链路在Matlab里写出来不超过200行。这篇文章就把这套方案的原理、代码、标定门限的坑一次讲清楚适合正在做协作频谱感知课程设计、毕业设计或者通信仿真的读者。1. 先搞懂协作频谱感知里数据到底怎么“融合”1.1 单节点检测的痛点和集中式融合的必然性我记得刚接触频谱感知那会儿最直接的想法就是每个节点做能量检测算一下接收信号的平均功率跟预设门限比一比大于门限就认为主用户存在。看起来很简单可真放到实际场景里单节点检测的问题特别多。最典型的是隐藏终端问题——次用户被高大建筑物挡住接收信号衰减很厉害明明主用户正在发射本地节点却测不出能量。其次是阴影衰落和频率选择性衰落会让单个节点的信噪比瞬间掉下去检测概率跟着崩。这时候就需要协作频谱感知。多个次用户分布在空间不同位置它们同时经历深度衰落的概率要比单个节点小得多。把这些节点的信息汇总起来再判决相当于利用了空间分集增益能明显改善检测性能。协作频谱感知按融合方式可以分为集中式、分布式和中继式其中集中式最常用、也最容易在Matlab里仿真每个次用户SU把本地感知结果通过控制信道发给融合中心FC融合中心对所有结果做数据处理最终输出一个全局判决。集中式融合又分硬融合和软融合。硬融合是每个节点先做本地判决上报1bit或者2bit的判决结果融合中心用AND、OR或者K-N规则汇总软融合则是节点把原始统计量比如能量值、匹配滤波输出直接发给融合中心融合中心用等增益合并、最大比合并等方式处理。硬融合通信开销小但本地判决已经丢了信息硬融合再怎么优化也有天花板软融合理论上保留更多信息性能上限更高。我们今天要说的PR指数检测器就属于软融合而且不是简单的加权平均是一种基于分布距离的融合策略。1.2 融合中心的判决本质这是在做一个分布检测很多教材在讲协作感知时会把融合中心简化为“把收到的统计量加起来再跟门限比”。这个思路没有错但它默认了一个前提所有节点的统计量服从同一个已知分布而且合并后的统计量仍有简单的理论分布。现实里这个前提往往不成立。不同节点所在位置的噪声功率不一样信道衰落系数不一样甚至采样点数也可能不一样直接把能量值加起来并不稳健。换个角度看融合中心拿到的其实是一堆来自不同节点的统计量样本。假如主用户不存在这些统计量应该符合“只有噪声”时的分布假如主用户存在一部分节点收到的信号会被主用户信号污染统计量的整体分布就会偏离纯噪声分布。于是判决问题变成了一个分布检测问题我手头这组统计数据和“纯噪声基准分布”差多少差得足够远就认为有主用户信号。这种思路的好处是不依赖所有节点完全同分布也不依赖每个节点的信噪比信息。只要主用户信号让足够多的节点统计量发生变化融合中心就能通过“分布偏离程度”感知到异常。我们下面要说的Pietra-Ricci指数检测器正是用来量化这个“偏离程度”的。2. Pietra-Ricci指数检测器的核心细节与选型逻辑2.1 指数到底是什么从距离定义说起Pietra-Ricci指数在信号处理领域并不算大热门但它有一个非常漂亮的定义对于两个概率密度函数p(x)和q(x)PR指数定义为二者总变差距离Total Variation Distance的二分之一也就是[ PR(p,q)\frac{1}{2}\int_{-\infty}^{\infty} |p(x)-q(x)| dx ]这个式子的几何意义很直观把两条概率密度曲线画在一个坐标系里它们重叠区域越大积分值越小完全不重叠时积分值为1完全相同时积分值为0。所以PR指数的取值被限制在[0,1]之间天然适合做检测门限。有人可能会问为什么要用PR指数而不是KL散度或Bhattacharyya距离我在对比过几种度量之后觉得PR指数有两个工程上很舒服的特性。第一它是对称的(PR(p,q)PR(q,p))而KL散度不对称分析时容易踩坑。第二它对分布尾部不敏感即使经验分布里出现个别极端样本积分也不会爆炸相比之下KL散度在两个分布支持下不相交时会趋向无穷大在统计量样本较少时非常不稳定。Bhattacharyya距离虽然也是对称且有界但它需要对密度函数做平方根运算用核密度估计时对带宽参数更敏感稍不注意就出现数值抖动。PR指数只做绝对值和积分简单、稳定、可解释性强所以我更推荐在频谱感知融合里用PR指数。2.2 用于频谱感知时统计量怎么构造用PR指数做集中式融合第一步是确定每个节点上报给融合中心的统计量。最常见的选择是能量统计量。假设第i个节点的接收信号为(y_i(n))采样点数为N那么本地能量统计量可以写成[ E_i\frac{1}{N}\sum_{n1}^{N}|y_i(n)|^2 ]这里的(y_i(n))通常是复基带信号。如果只有复高斯噪声且噪声方差为(\sigma^2)那么(|y_i(n)|^2)的均值为(\sigma^2)方差也是(\sigma^2)。按大数定律当N足够大时(E_i)近似服从均值为(\sigma^2)、方差为(\sigma^4/N)的高斯分布。为了方便后续统一建模可以对能量做归一化[ T_i\frac{E_i}{\sigma^2} ]这样在纯噪声假设H0下所有节点的(T_i)都近似服从均值为1、方差为(1/N)的高斯分布。这个归一化过程在Matlab里特别重要很多仿真结果对不上十有八九是忘了做这一层处理。融合中心收到K个节点的(T_1,T_2,\dots,T_K)之后不需要假设它们服从什么特定分布直接用这些样本得到经验累积分布函数ECDF。同时根据H0假设下的理论模型可以画出一条理论噪声CDF曲线即均值为1、标准差为(1/\sqrt{N})的高斯CDF。接下来计算这两条CDF之间的PR指数[ D\frac{1}{2}\int_{-\infty}^{\infty} \left|\hat{F}(x)-F_0(x)\right| dx ]其中(\hat{F}(x))是经验CDF(F_0(x))是理论噪声CDF。如果D大于门限(\lambda)融合中心判决主用户存在否则判决不存在。这里要特别说明一点为什么用CDF而不是PDF因为PDF需要做核密度估计带宽选择不当会直接毁掉结果而经验CDF直接由样本排序计算不需要调参收敛性也有保证。用数值积分计算两条台阶状曲线之间的面积Matlab的trapz函数就能轻松完成效率很高。2.3 为什么PR检测器能提升协作感知性能传统软融合方案如等增益合并或最大比合并本质上是把多个统计量加权平均。这种方案在理想条件下性能不错但有一个隐含假设参与融合的节点质量相近、没有特别离谱的异常值。实际环境里部分节点可能处于深度衰落也可能受到本地干扰上报的能量统计量会严重偏离正常范围。此时简单平均会被异常值带偏融合中心反而容易做出错误判决。PR指数检测器的思路完全不同。它考察的是整体分布偏离基准的程度而不是某个统计量的大小。即便一两个节点给出异常值它们对经验CDF的影响也有限不会像加权平均那样被直接“拉走”。换句话说PR指数对离群值更稳健。另外主用户信号出现时不只是均值变化统计量的整个分布形状都会变化PR指数能够捕捉这种形状变化而传统能量检测只盯着均值是否抬高丢掉了一些信息。还有一个实际好处PR指数检测器不需要知道各节点的信道增益或信噪比也不需要给每个节点分配权重。融合中心只需要一个“纯噪声统计量分布”的参考模型就可以完成判决。这在节点数量动态变化的场景里特别友好——节点增多时经验CDF会变得更平滑PR指数的积分结果更稳定节点减少时也不会出现算法崩溃。3. Matlab实现从场景建模到画出ROC曲线3.1 仿真场景和参数设置我习惯先定下一组基础仿真参数把主用户信号、信道、噪声和节点数都固定下来快速跑通整个流程再逐步调整参数看性能变化。下面的参数设置可以作为参考参数取值说明协作节点数K8可调为4/16节点越多融合效果越好每节点采样点数N1000能量估计精度相关主用户信号调制方式QPSK随机符号功率归一格信道类型瑞利平坦衰落每个节点的衰落系数独立噪声分布复高斯白噪声实部虚部方差各0.5信噪比范围-15dB到0dB关注低信噪比场景蒙特卡洛次数10000保证检测概率曲线平滑目标虚警概率0.1用于门限标定主用户信号我一般用QPSK因为它实现简单、带宽利用率高而且对能量检测来说具体的调制方式并不重要只要信号功率存在即可。每个节点经过独立的瑞利衰落信道再加上复高斯白噪声。这里要注意噪声功率的归一化设定如果信号平均功率设为1那么给定信噪比SNR_dB后噪声方差就是(\sigma_n^210^{-SNR_dB/10})。但为了后面PR指数计算简单我会让噪声方差固定为(\sigma_n^21)通过调整信号幅度来控制信噪比这样理论噪声CDF的均值和方差不用动态修改。3.2 核心代码实现能量统计量、经验CDF与PR指数计算先写一个计算PR指数的函数输入是K个节点的能量统计量T以及理论噪声分布的均值mu0和标准差sigma0输出是PR指数D。这里我用经验CDF和理论CDF的绝对面积差来近似积分。function D pr_index_from_energy(T, mu0, sigma0) % T : Kx1 向量每个节点上报的归一化能量统计量 % mu0 : H0假设下统计量的均值 % sigma0: H0假设下统计量的标准差 % D : Pietra-Ricci指数 % 经验CDF [f_emp, x_emp] ecdf(T); % 生成理论CDF需要覆盖经验CDF的有效区间 x_min min(x_emp); x_max max(x_emp); x_grid linspace(x_min, x_max, 500); % 对经验CDF插值到统一网格 f_emp_interp interp1(x_emp, f_emp, x_grid, previous, extrap); % 理论噪声CDF f_theory normcdf(x_grid, mu0, sigma0); % PR指数 0.5 * sum(|经验CDF - 理论CDF|) * 步长 D 0.5 * trapz(x_grid, abs(f_emp_interp - f_theory)); end这里用interp1的’previous’选项是刻意保留经验CDF的阶梯形状如果用线性插值会把台阶抹平PR指数会有偏差。实际上更好的做法是直接把经验CDF的每个跳跃点都保留下来和理论CDF在同样的采样点上比较但在大部分场景下上述近似已经足够。下面这段代码演示了怎么在H0假设下标定门限并在H1假设下计算检测概率。clear; clc; rng(1); % 参数设置 K 8; % 协作节点数 N 1000; % 每节点采样点数 Pfa_target 0.1; % 目标虚警概率 nMc 5000; % 蒙特卡洛次数 SNR_dB -10; % 信噪比 % H0下理论统计量分布参数 mu0 1; sigma0 sqrt(1/N); % 预分配 pr_null zeros(nMc, 1); pr_signal zeros(nMc, 1); % 第一步H0下估计PR指数的分布并标定门限 for mc 1:nMc noise (randn(K, N) 1i*randn(K, N)) / sqrt(2); E mean(abs(noise).^2, 2); % 每个节点能量 T E; % 噪声方差为1已经归一化 pr_null(mc) pr_index_from_energy(T, mu0, sigma0); end threshold quantile(pr_null, 1 - Pfa_target); % 第二步H1下计算检测概率 signal_power 10^(SNR_dB/10); detect_count 0; for mc 1:nMc signal (randi([0 3], 1, N) * 2 - 3) ... 1i*(randi([0 3], 1, N) * 2 - 3); signal signal / sqrt(mean(abs(signal).^2)); % 信号功率归一到1 signal sqrt(signal_power) * signal; % 按信噪比缩放 h (randn(1, K) 1i*randn(1, K)) / sqrt(2); % 瑞利衰落 % 每个节点接收信号 信号*衰落 噪声 Y zeros(K, N); for k 1:K Y(k, :) h(k) * signal (randn(1, N) 1i*randn(1, N)) / sqrt(2); end E mean(abs(Y).^2, 2); pr_now pr_index_from_energy(E, mu0, sigma0); if pr_now threshold detect_count detect_count 1; end end Pd detect_count / nMc; fprintf(SNR %.1f dB, Pfa %.3f, Pd %.3f\n, SNR_dB, Pfa_target, Pd);这段代码简洁但完整。第一层蒙特卡洛循环只产生纯噪声每次得到一个PR指数累计nMc次后用quantile函数取1-Pfa分位数作为门限。第二层循环在H1条件下计算检测概率。我在实际跑仿真时会把两层循环包成一个大函数并用parfor并行加速。3.3 门限怎么定蒙特卡洛标定不能省PR指数检测器的一个坑是很难直接推导出门限的闭式表达式。虽然理论上当K和N都趋近无穷时经验CDF会收敛到理论CDFPR指数趋近0但有限样本下的PR指数分布并没有简单的解析形式。所以最稳妥的方法就是蒙特卡洛标定——在H0假设下反复仿真得到PR指数在没有主用户时的经验分布取某个高分位数当门限。我在实验中观察到PR指数的分布跟K和N都有关系。K越大经验CDF越稳定PR指数整体越小分布也越集中N越大理论噪声CDF越精确同样会让PR指数变小。所以门限不是常数必须针对你实际的K和N重新标定。换一组参数就沿用旧门限虚警率大概率会翻车。标定门限的代码我习惯写成函数function thr calibrate_threshold(K, N, Pfa, nCal) pr_null zeros(nCal, 1); for mc 1:nCal noise (randn(K, N) 1i*randn(K, N)) / sqrt(2); E mean(abs(noise).^2, 2); pr_null(mc) pr_index_from_energy(E, 1, sqrt(1/N)); end thr quantile(pr_null, 1 - Pfa); endnCal建议至少取5000否则门限本身抖动也很大。如果仿真时间紧张可以先跑2000次粗标再逐步加密。我自己的经验是门限标定的蒙特卡洛次数不要低于检测概率仿真次数否则最终ROC曲线在高Pfa区域会出现不光滑的阶梯。3.4 画图展示结果从单点数字到完整ROC曲线单看一个SNR下的Pd根本看不出算法优劣。我建议至少画两种图ROC曲线Pd与Pfa的关系和Pd随SNR变化曲线。以总虚警概率0.01、0.05、0.1等一组取值为横坐标通过蒙特卡洛标定得到对应门限再在各个门限下计算Pd就能画出ROC曲线。代码如下Pfa_range 0.01:0.01:0.30; Pd_range zeros(size(Pfa_range)); for idx 1:length(Pfa_range) thr calibrate_threshold(K, N, Pfa_range(idx), 2000); % 计算H1下的PR指数并统计检测概率 pd_tmp 0; for mc 1:nMc % 生成信号、衰落、噪声计算PR指数与thr比较 % ... end Pd_range(idx) pd_tmp / nMc; end figure; plot(Pfa_range, Pd_range, b-o, LineWidth, 1.5); xlabel(虚警概率 Pfa); ylabel(检测概率 Pd); title(PR指数集中式融合的ROC曲线); grid on;ROC曲线越靠近左上角说明系统性能越好。我实测下来K从4增加到8同样SNR下检测概率大约能提升10到15个百分点再增加到16提升幅度变缓。这说明协作感知的增益并不是线性增长的设计系统时需要权衡节点数量和反馈开销。4. 常见问题与排查技巧实录4.1 为什么经验CDF和理论CDF的数值差很小导致检测概率低这是我第一次跑PR指数融合时遇到的第一个问题明明主用户信号已经很强了融合中心算出来的PR指数还是很小门限设得低一点虚警又压不住。后来排查发现问题出在能量统计量没有做归一化。我直接把每个节点的均值能量E_i送进计算函数H0下E_i的均值等于噪声方差(\sigma_n^2)方差是 (\sigma_n^4/N)。如果噪声方差不是1理论CDF还傻乎乎地按均值为1、标准差为1/N去算两条曲线自然对不上。解决方法是把能量统计量除以噪声方差让H0下的均值固定为1。如果噪声方差未知可以先在系统空闲时估计一段噪声功率再拿这个估计值归一化。即使估计有一点误差只要不是差得非常离谱PR指数的性能退化也不明显这是它比纯能量门限检测稳的一个重要原因。4.2 PR指数对门限太敏感虚警率不可控怎么办门限标定的核心是“和算法使用场景同分布”。如果你在标定门限时没有经过瑞利衰落只有纯高斯噪声而在实际检测时信号路径有衰落那门限本身问题不大因为门限只需要描述H0行为但如果你的接收机噪底变了或者采样点数变了门限就必须重新标定。还有一种情况是你用了ksdensity估计PDF然后计算PR指数。这时门限对核函数带宽极其敏感带宽太大分布被抹平成一条胖曲线PR指数偏小带宽太小又会出现一堆毛刺PR指数偏大。解决建议是尽量用经验CDF计算别碰核密度估计。CDF不需要调参退一万步讲就算出现异常样本台阶状CDF也只是多了一个小台阶不会像密度估计那样出现一个尖峰。如果发现虚警率还是压不住先画一下H0下PR指数的直方图看看是不是双峰或长尾。如果是大概率是随机数生成出了问题比如矩阵维度不对或者能量计算把复数模平方写漏了。4.3 仿真很慢、结果抖动大怎么办仿真慢通常是因为蒙特卡洛次数太大并且第二层循环里又套了一层遍历节点的for。我写的示例为了可读性保留了节点循环但实际工程中完全可以用矩阵运算一次性生成所有节点的数据。把噪声一次性生成K×N矩阵衰落系数生成1×K向量利用广播机制得到接收矩阵代码更简洁速度也更快。如果机器有多核强烈建议用parfor替换外层蒙特卡洛循环。我自己的经验是8核机器上并行大约能提速4-6倍10000次蒙特卡洛之前要跑一分钟并行后十几秒就能结束。要注意的是并行循环内不能再调用calibrate_threshold这种每次都跑几千次的函数否则每轮都在重复标定反而拖慢速度。正确做法是先把门限标定好再进入并行检测步骤。4.4 本地节点数量不等时经验CDF怎么处理实际系统中不同节点可能上报的统计量样本数不同比如有的节点采样点少能量估计方差大。经验CDF对这种异方差问题并不敏感因为PR指数只关心整体分布的形状差异不关心单个样本的方差。不过如果某个节点采样点数特别少它的能量统计量会出现很大抖动经验CDF尾部会拉出很长的尾巴导致PR指数虚高。遇到这种情况我建议给每个节点设置最小采样点数或者用加权经验CDF让高可靠性节点贡献更大的权重。这个优化方向在论文里可以作为延伸工作但在平时课程设计和工程验证中只要保持所有节点采样点数一致就能避免大部分问题。5. 扩展思路这套方法还能怎么用PR指数检测器并不局限于能量统计量。实际上任何能区分H0和H1的本地特征量——比如循环平稳特征、协方差矩阵特征值、匹配滤波输出——都可以作为节点上报的统计量只要H0下的理论分布可以近似得到。融合中心照样能构造经验CDF和理论CDF的比较计算PR指数。我试过用信号循环谱的峰值作为本地统计量在cyclostationary特征检测场景下PR指数融合的抗噪声性能比能量检测更高代价是计算量变大。如果你正在做的项目不要求实时性这个方法完全可以当成一个对比算法。另外如果融合中心到节点之间的控制信道会出错上报的统计量可能丢包或被篡改。经验CDF对少量丢包并不敏感因为缺失几个样本不会显著改变分布形状。这一点比硬判决融合更有优势——硬判决丢一个关键节点的1bit可能就影响最终结果。最后分享一个小习惯每次跑完PR指数检测器仿真的结果我会把H0分布、H1分布和门限三条线画在同一张图上。这张图不仅是论文里的素材也是自查逻辑的利器。只要看见门限位置明显不合理比如落在H0分布峰值上就知道参数设置肯定有问题。仿真调试这件事光靠看最终数字远远不够眼见为实永远是最快的排查方式。