最近在做认知无线电的频谱感知仿真时又翻出了Pietra-Ricci指数检测器。这个检测器在论文里出现的频率不如能量检测器高但它在集中式协作感知里真的很能打。尤其是当你面临噪声功率不确定、各节点信噪比参差不齐的场景时用Pietra-Ricci指数来做数据融合比单纯把能量值相加或者用硬判决合并的思路要稳得多。这篇内容围绕“Pietra-Ricci指数检测器用于集中式数据融合协作频谱感知”展开我会把原理、系统模型、Matlab实现、仿真细节以及我踩过的坑全部梳理一遍。无论你是准备课程设计还是想给自己的论文补充仿真结果这份实操记录都能直接参考。1. 为什么要用协作频谱感知单节点检测的局限与融合的价值1.1 频谱感知解决什么问题认知无线电的核心任务之一就是让次用户在不干扰主用户的前提下机会式地接入空闲频谱。这就需要一个可靠的“先听再说”机制也就是频谱感知。常见做法是次用户在某个频段上接收一段信号通过能量检测、匹配滤波检测或者循环平稳特征检测等方法判断主用户信号是否存在。判断结果可以写成二元假设检验H0主用户不存在接收信号只有噪声。H1主用户存在接收信号是主用户信号经信道传播后的结果加上噪声。这个二元判决看似简单但实际场景里单个次用户的感知性能非常脆弱。一旦出现阴影衰落、多径衰落或者突发干扰本地检测结果可能偏离真实情况导致虚警或者漏检。1.2 单用户检测为什么靠不住单个节点做能量检测时需要知道噪声功率的准确值。如果噪声功率估计偏差了3dB检测性能可能直接崩掉。更麻烦的是单个节点可能正好处于主用户信号的深衰落位置接收信号强度极低无论怎么提高检测门限都救不回来。在真实环境里主用户到各个次用户的传播路径不同有的次用户离得近接收信噪比高有的次用户被建筑物挡住几乎收不到信号。如果只依赖某一个用户的判断结果就完全赌在这个节点的信道上性能很容易出现较大波动。协作频谱感知的好处是把多个次用户的信息汇总起来利用空间分集来对抗衰落。即使部分节点处于深衰落只要其他节点的信号质量尚可整体判决仍然可靠。这也是认知无线电研究中协作感知被广泛采用的根本原因。1.3 集中式数据融合与分布式融合的取舍协作感知的数据融合方式主要分两种集中式融合各次用户把本地观测数据或本地统计量发送到融合中心融合中心统一做出全局判决。分布式融合各节点先本地判决再把1bit结果发送给邻居或融合中心最后按某种投票或者加权规则合成。集中式的优点是信息保留完整理论上性能优于分布式。缺点是传输开销大并且依赖融合中心。分布式虽然开销小但本地硬判决会丢失大量信息性能增益有限。本文讨论的Pietra-Ricci指数检测器属于集中式融合的一种特殊实现它不需要各节点上报原始采样数据只需要各节点上报本地能量值。这在带宽受限的认知无线电网络里非常讨巧既保留了软信息又避免了巨大的传输开销。2. Pietra-Ricci指数检测器原理拆解2.1 从“统计差异”到“信号存在性”Pietra-Ricci指数最早并不是为通信设计的。它来自经济学和统计学用来衡量收入分配的不平等程度。想象一个社会有K个人每个人的收入为(x_i)。如果所有收入完全相等那么收入分布就是绝对平均的如果一部分人占有了几乎所有收入收入分布就高度不平均。Pietra-Ricci指数正是量化这种不平均的指标之一。把它搬到频谱感知里K个次用户的本地能量值(e_1, e_2, ..., e_K)就相当于K个人的收入。当主用户不存在时所有次用户接收到的都是纯噪声。尽管噪声是随机的但各用户能量统计量围绕同一个噪声功率水平波动整体上差异不大PR指数会很小。当主用户存在时主用户信号经过不同路径到达各次用户由于路径损耗、阴影衰落的存在各用户接收到的信号功率明显不同。能量值之间的差异被拉大PR指数也随之变大。所以用PR指数做判决统计量本质上是在问这一组能量值的分布是否足够“不均匀”如果足够不均匀就倾向于认为背后存在一个共同的主用户信号源只是因为信道差异让信号功率分布散开了。如果能量值都很均匀那大概率大家收到的都是噪声。这种思路最大的好处是不需要知道噪声功率的绝对值。因为PR指数是一个归一化的相对差异度量对噪声功率的尺度变化不敏感。这对实际系统来说意义很大。2.2 PR指数的数学定义与计算细节对于一组非负数据(x_1, x_2, ..., x_K)Pietra-Ricci指数的离散估计可以写为[ PR \frac{1}{2K\bar{x}}\sum_{i1}^{K}|x_i - \bar{x}| ]其中(\bar{x} \frac{1}{K}\sum_{i1}^{K}x_i)是样本均值。这个式子的直观含义是把每个值与均值的绝对偏差累加起来再除以两倍的总和。它度量的是数据围绕均值的“平均绝对离差”占均值本身的比例。当所有(x_i)相等时PR0当数据极不均衡时PR趋近于1。这个公式与基尼系数有关但计算比基尼系数简单。它不涉及排序后的累加权重直接用绝对偏差就能算非常适合在Matlab里用向量化操作实现。在实际计算时需要留意几个细节如果某个用户的能量值因为数值下溢变成0公式不会出问题因为均值通常不会为0。如果K比较小比如只有2或3PR指数的取值分辨率很有限。举个例子K2时PR指数要么是0要么是某个非零固定值很难做精细的门限控制。因此实际使用中K不宜太小建议至少4个以上次用户参与融合。PR指数计算的是能量值之间的相对差异不关心能量值的绝对大小因此不同噪声功率水平下的用户组可以放在一起比较。2.3 检测统计量与判决规则整个检测流程可以写成每个次用户对接收信号(y_i(n))做能量积累得到本地能量值(e_i \frac{1}{N}\sum_{n1}^{N}|y_i(n)|^2)。各次用户将(e_i)上报给融合中心。融合中心收集到向量(\mathbf{e} [e_1, e_2, ..., e_K])。计算PR指数作为全局检验统计量。将PR指数与门限(\lambda)比较若(PR \lambda)判为H1否则判为H0。门限(\lambda)的设定通常基于奈曼-皮尔逊准则在给定虚警概率(P_{fa})约束下寻找令检测概率最大化的门限。由于PR统计量的解析分布很难求工程上最常用的是蒙特卡洛仿真法——在H0假设下多次仿真得到PR统计量的经验分布再取对应的分位数作为门限。3. 系统模型与仿真参数设计3.1 信号模型与信道假设仿真时我采用如下参数主用户信号使用BPSK调制符号序列随机生成信号功率固定为(P_s)。次用户数(K)一般取4、8、12。每个次用户的感知采样点数(N)例如1000。信道瑞利平坦衰落每个用户的主用户信道增益为复高斯随机变量归一化功率为1。也就是说(h_i \sim \mathcal{CN}(0,1))表示主用户信号到达第i个用户时经历了不同的衰落。噪声复加性高斯白噪声均值为0方差(\sigma_n^2)仿真时统一设噪声功率为1通过调整信号功率来控制信噪比。在H0情况下第i个用户接收信号为[ y_i(n) w_i(n) ]在H1情况下[ y_i(n) h_i \cdot s(n) w_i(n) ]其中(s(n))为主用户信号(w_i(n))为噪声。这里信道只考虑慢衰落即在一个感知周期内(h_i)保持不变。实际系统里如果快衰落使信道在感知窗口内快速变化能量值会有所平滑但PR指数的基本行为不变。3.2 融合中心如何利用本地能量值各用户上报的能量值(e_i)本质上是一个非负实数。融合中心拿到这K个数后直接套公式算PR指数。这里有几个值得注意的地方如果各用户上报的是平均能量即除以(N)后的能量PR指数不受(N)影响因为均值会消掉(N)的倍数。但如果各用户上报的是累加能量那么累加点数(N)会同时放大分子分母PR指数的值并不变。所以无论是否除以N计算出的PR指数都一致。实际工程中上报的软信息需要量化量化误差会影响PR指数。仿真中可以暂时忽略量化但如果你要做更贴近实际系统的性能评估建议给能量值加上均匀量化噪声观察性能损失。融合中心不要求各用户的噪声功率相等。PR指数的优势在于它比较的是能量分布的不均匀性而不是能量绝对大小所以不同用户的噪声系数不同也能工作。这一点在仿真中可以直接验证把不同用户的噪声方差设为不同值PR检测器的性能衰减少于传统能量检测。3.3 门限获取蒙特卡洛仿真与理论近似PR统计量的理论分布较难推导工程中直接通过蒙特卡洛仿真获取门限。做法是在H0条件下独立生成(M)次仿真建议(M10000)次以上。每次计算得到一个PR值共得到(M)个PR样本。将所有PR样本从小到大排序。给定目标虚警概率(P_{fa})取排序后索引为(\lceil M(1-P_{fa})\rceil)的样本值作为门限(\lambda)。因为PR指数存在上限大约不超过1且H0条件下的PR值集中在较小数值范围这个经验门限的稳定性主要取决于(M)。我的经验是(M)至少10000否则ROC曲线尾巴抖动很大。如果你只需要一个固定工作点可以取5000次但为了画ROC曲线建议多跑一些。另外也可以计算理论门限的近似值比如利用皮尔逊分布拟合H0下的PR统计量但实际价值不大。蒙特卡洛法在离线仿真中已经够用了。4. Matlab代码实现从零搭建仿真链路4.1 主程序框架与参数设置先给出主程序的结构再逐段解释。整个仿真分四块参数设置、门限估计、信号生成与判决、性能统计和绘图。参数设置代码%% 参数设置 clear; close all; clc; K 6; % 次用户数量 N 1000; % 每个用户感知采样点数 M 10000; % 蒙特卡洛次数 SNR_dB -15:2:5; % 信噪比范围dB Pfa 0.1; % 目标虚警概率 % 固定随机种子保证可复现 rng(2025);这里的K和N对检测性能影响很大。K越大PR指数的统计分辨率越高融合多样性增益越大N越大每个用户的能量估计越准PR值受到的不确定度越小。4.2 H0条件下门限估计函数门限估计需要在H0条件下生成能量矩阵。这里我们不依赖真实噪声功率只假设噪声功率为1。为了方便扩展我把噪声功率设为全局变量。sigma2 1; lambda zeros(1, length(SNR_dB)); % 每个SNR都可以有自己的门限因为门限在不同条件下会变但实际系统不知道SNR所以通常应该用一个统一门限这里需要小心如果门限随着SNR变化那就相当于开挂实际系统中SNR未知。正确做法是假设一个标称噪声功率按H0仿真得到门限然后在各种SNR下都用这个门限。由于我们设噪声功率固定为1信号功率变化只是改变SNRH0情况下没有信号所以H0统计量分布与SNR无关门限也应该相同。因此门限只需要算一次。% 在H0条件下估计门限 pr_samples zeros(M, 1); for m 1:M Y sqrt(sigma2/2) * (randn(K, N) 1i*randn(K, N)); e mean(abs(Y).^2, 2); % 每个用户的平均能量 pr_samples(m) PR_index(e); end pr_samples sort(pr_samples); lambda pr_samples(ceil(M * (1 - Pfa)));需要注意的是这里生成了(K \times N)的复噪声矩阵每个用户的平均能量e是一个K x 1的列向量。4.3 PR指数计算函数PR指数的函数实现非常简单function pr PR_index(e) % e: K x 1 非负能量向量 mu mean(e); pr sum(abs(e - mu)) / (2 * K * mu); end这个函数可以直接保存为PR_index.m。计算时注意分母中的K不能漏掉。4.4 信号生成、融合判决与性能统计主程序中对每个SNR进行蒙特卡洛仿真统计检测概率。Pd zeros(length(SNR_dB), 1); for s 1:length(SNR_dB) snr 10^(SNR_dB(s) / 10); Ps snr * sigma2; % 信号功率保证接收端平均信噪比为snr detected 0; for m 1:M % 生成主用户信号BPSK符号 s_tx 2 * randi([0 1], 1, N) - 1; % BPSK幅度为±1 s_tx sqrt(Ps) * s_tx; % 调整信号功率 % 每条用户信道 h (randn(K, 1) 1i * randn(K, 1)) / sqrt(2); noise sqrt(sigma2/2) * (randn(K, N) 1i*randn(K, N)); % H1: 主用户信号存在 Y h .* s_tx noise; % 本地能量检测 e mean(abs(Y).^2, 2); % 融合中心计算PR指数 pr PR_index(e); % 判决 if pr lambda detected detected 1; end end Pd(s) detected / M; end这段代码的运算量比较大因为H0门限估计和每个SNR的H1仿真都要跑上千次。实际运行时M和SNR点数可以适当调小比如M5000SNR范围为-10:2:0这样就比较快。4.5 ROC曲线绘制ROC曲线可以用不同的门限遍历来画。做法是在H0和H1条件下分别得到两组PR统计量样本然后对一组门限值计算对应的虚警概率和检测概率。% 生成H0和H1条件下的PR统计量以某个SNR为例 M 10000; pr_h0 zeros(M, 1); pr_h1 zeros(M, 1); snr -5; % 示例信噪比 Ps 10^(snr/10) * sigma2; for m 1:M % H0 Y0 sqrt(sigma2/2) * (randn(K, N) 1i*randn(K, N)); pr_h0(m) PR_index(mean(abs(Y0).^2, 2)); % H1 s_tx (2 * randi([0 1], 1, N) - 1) * sqrt(Ps); h (randn(K, 1) 1i * randn(K, 1)) / sqrt(2); noise sqrt(sigma2/2) * (randn(K, N) 1i*randn(K, N)); Y1 h .* s_tx noise; pr_h1(m) PR_index(mean(abs(Y1).^2, 2)); end % 门限遍历 thresholds linspace(0, 0.2, 200); Pfa_curve zeros(size(thresholds)); Pd_curve zeros(size(thresholds)); for t 1:length(thresholds) Pfa_curve(t) mean(pr_h0 thresholds(t)); Pd_curve(t) mean(pr_h1 thresholds(t)); end figure; plot(Pfa_curve, Pd_curve, LineWidth, 1.5); xlabel(虚警概率 P_{fa}); ylabel(检测概率 P_{d}); title([PR指数检测器 ROC曲线, SNR , num2str(snr), dB, K , num2str(K)]); grid on;绘制出的ROC曲线通常随着SNR提高向左上角移动。重点不是曲线本身多漂亮而是理解不同参数下曲线的变化趋势。5. 关键结果与性能对比5.1 不同信噪比下的检测概率变化在K6、N1000、Pfa0.1的条件下可以看到当SNR从-15dB增加到-5dB的过程中检测概率逐渐上升并趋近于1。PR检测器在SNR低于-10dB时性能下降得比较快但相比传统的单节点能量检测它的优势在于节点间的空间分集补偿了低SNR的劣势。需要提醒的是PR检测器并不是在所有场景下都优于能量检测。如果所有次用户经历完全相同的信道比如所有用户都紧挨着且没有衰落那么有信号时各用户能量也几乎相同PR指数并不会明显变大检测性能会很差。这也是PR检测器的适用边界它依赖节点间信道差异。5.2 用户数量对融合检测的影响用户数K增加时PR指数的计算样本增多统计分布更稳定门限也更容易区分H0和H1。但K增大到一定值后检测性能提升速度变缓。原因在于PR指数反映的是能量分布的不均匀程度当K很大时个别差异对整体指标的贡献被稀释所以融合增益存在饱和效应。仿真时建议分别跑K4、K8、K12对比曲线后你就能看到这种饱和趋势。选K时也要考虑实际系统的开销K越大上报链路越多网络流量和时延都增加。理性选择是8个左右。5.3 与能量检测和MME方法对比传统的集中式能量检测融合规则是把所有用户的能量求平均然后与噪声功率相关门限比较。这个方法在噪声功率准确已知时性能不错但噪声不确定性会带来平台效应也就是无论SNR怎么提高检测性能都无法超过某个上限。PR检测器天然免疫噪声功率水平的变化因为PR指数是归一化相对量。另一个常见的融合检测是基于样本协方差矩阵的最大最小特征值之比MME。MME同样不依赖噪声功率它的理论基础来自随机矩阵理论性能在相关信号场景下很好。但MME需要上传原始采样数据或者至少需要融合中心可以获得高维协方差矩阵通信开销大。PR检测器只上传一个标量计算复杂度也低很多。总结下来PR检测器的定位是在不要求极高检测精度的场景下用一个计算简单、对噪声不敏感、传输开销小的方案实现稳健的协作感知。6. 踩坑记录与实用建议6.1 蒙特卡洛仿真的收敛性与次数设置蒙特卡洛次数M如果不够PR统计量的经验分位数波动很大导致虚警概率偏离目标。我第一次用2000次结果在同一参数下重复跑几次门限值明显不同ROC曲线也不光滑。后来把M提到10000以上才稳定下来。另外门限估计和性能统计要使用不同的随机序列。如果都用同一组随机数据会产生一种“过拟合”效应虚警率看起来很好但换一组数据就露馅。建议在门限估计和性能仿真之间重新设置随机种子或者缓存两组独立数据。6.2 噪声不确定性的处理PR检测器对绝对噪声功率不敏感但如果不同用户的噪声功率本身就不同比如用户1的噪声方差是1用户2的噪声方差是2PR指数的行为会发生变化。仿真时可以故意设置不一致的噪声方差观察检测概率变化。我看到的结果是PR检测器依然能工作但门限需要按实际噪声统计来校准。如果你的场景里各用户噪声功率差异巨大建议先做噪声归一化比如让每个用户上报“归一化能量”。6.3 用户间相关性导致PR指数失真如果多个次用户距离很近它们接收到的信道衰落高度相关那么有信号时能量差异反而可能减小PR指数被压低。这在多天线节点或者密集部署场景里尤其明显。解决办法是尽量选择空间分布分散的用户参与协作或者在融合时给用户能量加权重但这又引入了额外的参数。仿真时可以用相关衰落信道模型来检验这种退化程度。6.4 门限偏移与动态调整在真实系统中噪声统计可能随时间缓慢变化。PR指数的门限如果固定不变实际虚警概率会偏离设计值。一种实用的校正是周期性地利用无信号时隙重新估计门限例如在认知无线电开窗检测的静默期内运行H0采样。Matlab仿真的话可以把门限估计写成一个函数模块方便替换成自适应门限策略。7. 扩展方向与个人体会PR检测器这个思路值得继续深挖的方向有不少。比如把PR指数和能量检测结合做双门限或者联合检测或者在分布式融合中用PR指数对节点可靠性进行加权剔除异常节点。我试过在融合中心对能量值先做一次排序再计算截断PR指数也就是只保留中间一半的值这样能抑制极端异常上报性能更稳。另外Pietra-Ricci指数的计算复杂度极低很适合部署到资源受限的嵌入式设备上。如果未来想以这种方式实现实时协作感知只需要把能量值通过控制信道传回融合中心存储和计算压力都很小。我在实际使用中发现这类非参数化的融合检测器最大的优点在于“不敏感”。它不像能量检测那样害怕噪声功率估计偏差也不像匹配滤波那样要求主用户信号结构完全已知。对于工程实现而言这是一种非常省心的选择。如果你正在做频谱感知的仿真或者课程设计建议不要只跑出几张图就结束。试着改变信道模型、用户数、采样点数观察PR指数统计量如何变化然后你再回头看这个检测器的设计思路会很自然地理解为什么它能在协作感知场景里站稳脚跟。代码本身并不复杂复杂的是搞懂每一个参数到底在改变什么。跑通之后把H0门限估计、H1性能统计、ROC曲线绘制拆成独立的测试模块后续再对比其他检测器就是水到渠成的事了。
