做协作频谱感知研究的朋友大概率都绕不开一个场景多个认知节点各自感知频段把统计量汇到融合中心做联合判决。这篇文章想聊的是基于Pietra-Ricci指数检测器的集中式数据融合协作频谱感知方案以及配套的Matlab实现。这种检测器利用样本协方差矩阵特征谱与理论噪声谱的偏离度来判定主用户是否存在天然不依赖噪声功率精确估计在多节点协作时又能通过数据融合换检测增益很适合正在做认知无线电、频谱感知算法对比或者课题复现的读者。我会把原理、融合策略、门限设计、仿真代码和踩过的坑一次讲完。1. Pietra-Ricci指数到底在度量什么1.1 一句话理解Pietra-Ricci指数Pietra-Ricci指数简称PR指数本质上是两个概率分布之间“偏离程度”的一种度量。它最早在经济学里用于衡量收入分布与均匀分布的差距后来被引入信号处理用来量化观测数据分布与假设分布的一致性。给定两个概率密度函数 p(x) 和 q(x)PR指数的定义是 ( P \frac{1}{2} \int |p(x)-q(x)| dx )。在离散情况下就是对应概率质量差的绝对值之和再乘 1/2。取值范围在0到1之间0表示两个分布完全重合1表示两个分布完全不重叠。不必纠结它的数学出处关键是抓住直觉这个指数统计的是“两个分布错开的部分”占总体积的比例。如果用一个生活化类比想象两个班学生的身高分布PR指数近似于回答“要从两个班里随机各抽一个人光看身高有多大概率能判断出他来自哪个班”。两个分布分得越开指数越大。这个特性放到频谱感知里尤其好用因为检测问题本质上就是区分“只有噪声”和“噪声加信号”两个状态状态之间分得越开判决就越可靠。1.2 从统计距离到频谱感知特征把PR指数用到频谱感知上核心思路是构造两个需要区分的分布。以多天线或过采样接收为例在纯噪声假设 H0 下接收信号样本协方差矩阵的特征值会聚集在一个相对确定的理论范围内。如果存在主用户信号 H1协方差矩阵的低秩分量会被“顶”出来特征值谱整体会偏离原来的噪声谱。PR指数检测器做的事情就是把“观测到的特征值分布”和“理论噪声特征值分布”之间的PR距离作为检测统计量。距离越大越倾向于判为有信号。之所以选择特征谱而不是直接用接收能量做统计量是因为能量检测很容易被噪声功率不确定性干扰。噪声功率一旦估计偏了能量检测的门限就可能失效而基于特征谱的PR指数检测器不需要把噪声功率作为先验输入只需要特征值之间的相对关系或者经验谱与理论谱的分布差异。这种特性在协作频谱感知里尤其有价值因为多个节点的噪声底噪不完全一致要逐个估计噪声功率反而麻烦。用PR指数每个节点都能独立输出一个“谱偏离程度”融合中心再对这些统计量做集中式处理形成全局判决。因此PR指数检测器的定位不是替代所有检测器而是在低信噪比、噪声不确定性明显、多节点可协作的场景下提供一种稳健的候选方案。2. 集中式协作频谱感知的系统设计与选型2.1 协作感知为什么必要单节点频谱感知的软肋在于信道效应的随机性。一个认知节点可能因为深衰落、建筑物遮挡或者多径抵消完全接收不到主用户信号出现“隐藏终端”问题。这时候无论本地检测器多聪明输入信号本身就是低的输出自然不可靠。协作感知指望的是空间分集多个节点分布在不同的地理位置同时陷入深度衰落的概率远低于单一节点。只要有一个节点感知到主用户信号融合中心就能从全局角度做出“信道被占用”的判断从而降低漏检概率。协作感知的代价是通信开销和融合复杂度增加。每个节点需要把本地感知结果或者原始数据上报给融合中心这就产生了控制信道的占用以及数据同步、衰落信道上报差错等问题。因此设计得当的协作感知需要在检测增益和开销之间做权衡。集中式数据融合结构可以理解成“中心决策树”模式所有节点是树叶融合中心是树干。节点只负责采集和初步处理决策权归中心。2.2 集中式数据融合的两种常见结构第一种结构是“先检测后融合”。每个节点用PR指数检测器算出一个标量统计量 (T_i)通过控制信道传给融合中心融合中心对 (T_i) 做等增益合并EGC或加权合并再与全局门限比较。这种结构最省开销节点只需上传一个浮点数而且对控制信道带宽要求很低适合节点数较多、信道资源紧张的分布式网络。第二种结构是“先上传后检测”。每个节点不再只上传统计量而是上传自己估计的样本协方差矩阵或者上传从协方差矩阵提取的特征向量与特征值。融合中心把所有节点信息拼成一个总体的观测矩阵再统一计算PR指数。这种结构保留了更多信息能够处理节点间信号相关的情况但上传量随天线数和采样数增加对同步要求也更高。我个人的建议是课题仿真阶段两种都写一下但重点放在第一种因为它在工程上更容易落地性能也已经足够说明协作增益。真正需要第二种结构的场景通常是节点本身有多根天线而且主用户信号在多个节点间相关这时候“先上传后检测”能获得联合处理的额外增益。标题里提到的“集中式数据融合”更多时候指的是融合中心把各节点上报的本地统计量汇总起来做全局判决所以我下面仿真也围绕第一种展开同时留一个第二种的扩展注释。2.3 信道模型与同步问题仿真建模时不能只看理想高斯信道。常见的做法是采用带路径损耗和块衰落的信道模型每个节点到主用户之间的信道系数 (h_i) 服从复高斯分布幅度遵循瑞利衰落节点噪声之间彼此独立。主用户信号如果是OFDM或者QPSK调制可以建模成复基带信号接收端采样N个点。PR指数检测器对同步误差其实并不算特别敏感因为它是基于协方差矩阵的累积统计量轻微的采样定时偏差只会让信号能量分散到更多特征值上不会像循环平稳检测那样直接丢失特征。不过如果采样偏差导致符号间干扰严重检测性能还是会有可感知的下降。因此在集中式融合的仿真里我建议把信道随机关心一次作为一个“慢变块”处理每个蒙特卡洛循环独立生成信道和噪声。同步因素可以先用理想同步简化随后在系统级仿真里再引入采样相位偏移观察PR指数统计量的变化范围。这样能分离算法增益和实现损耗。3. 检测统计量构造与门限怎么定3.1 从特征值到PR指数统计量这里需要给出具体的统计量构造过程。设第 i 个节点接收 (N_s) 个采样点、(M_a) 根天线或 (M_a) 个过采样支路接收矩阵为 (X_i)则样本协方差矩阵 (R_i X_i^H X_i / N_s)维度 (M_a \times M_a)。对 (R_i) 做特征值分解得到特征值 (\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_{M_a})。理论噪声特征值分布可以用Marchenko-Pastur律近似为了工程实现常用的简化处理是先把特征值除以均值归一化再将归一化特征值排序计算它们的经验累积分布函数 (F_{obs}(x))并与理论累积分布 (F_{noise}(x)) 做差。PR指数检测统计量可以写成( T_{\text{PR}} \frac{1}{2M_a} \sum_{k1}^{M_a} | F_{obs}(\lambda_k) - F_{noise}(\lambda_k) | )这里 (\lambda_k) 是归一化后的特征值(F_{obs}(\lambda_k) k/M_a)。实际写Matlab时可以用排序后的特征值直接构造经验CDF序列再和指数分布的理论CDF逐点求差。为什么可以和指数分布比因为当样本数足够大、天线数固定时单个样本协方差矩阵的特征值在纯噪声下会围绕总体均值波动归一化后近似服从某种集中在1附近的分布具体形状与 ((M_a, N_s)) 有关。严格做法是查表或仿真得到该分布仿真阶段用蒙特卡洛门限更稳妥。这里有几种常见检测器的对比方便看清楚PR指数检测器的定位检测器是否依赖噪声功率先验主要优势主要劣势能量检测依赖实现简单受噪声不确定性影响最大最小特征值检测不依赖稳健小样本下门限理论不完善PR指数检测器不依赖对分布整体偏离敏感协作友好需要计算特征值分解复杂度略高3.2 门限设计CFAR思想检测器性能最终靠门限体现。基于PR指数的统计量在H0下的分布没有标准闭式表达式尤其在有限采样点数下。我在实际复现时用的是“蒙特卡洛标定门限”的做法保持所有参数不变只生成纯噪声数据重复 (M_c) 次得到 (M_c) 个 (T_{\text{PR}}) 值排序后取 ((1-P_{fa})) 分位数作为判决门限。比如目标虚警概率0.01蒙特卡洛10000次就取排序后第9900个值。这样做虽然仿真时间会长一点但门限完全贴合实际分布避免理论近似带来的系统性偏差。要小心的一点是门限必须和统计量在同一条件下标定采样点数、天线数、节点数、融合方式变了H0分布都会变。有些同学把单节点门限直接拿去做四节点融合结果检测概率虚高因为融合后的统计量方差变小原来的门限偏大或偏小都不一定。正确做法是每个仿真参数组都独立标定或至少用足够多的噪声样本拟合分布。还可以在仿真里做一个“门限-虚警率校准曲线”确认实际虚警率接近目标值再测检测概率否则后面画ROC就是自欺欺人。4. Matlab仿真单节点检测与集中式融合全流程4.1 仿真场景参数配置仿真我按下面的参数跑通了整个流程你可以直接套用或修改参数取值说明主用户信号QPSK复基带信号频谱感知常用调制采样点数 (N_s)1000每节点每蒙特卡洛采样天线数 (M_a)2节点本地可配多天线节点数 (K)1, 2, 4, 8对比协作增益信噪比范围-20 dB ~ 0 dB重点看低信噪比区蒙特卡洛次数5000跑门限时用10000虚警概率目标0.01CFAR门限信道模型瑞利块衰落每MC独立生成注意Pietra-Ricci指数检测器本身可以工作在单天线节点但为了体现“特征谱”的优势至少让每个节点有2个接收支路可以是2根天线也可以是对同一信号做2倍过采样。如果节点只有单天线样本协方差矩阵退化成一个标量PR指数就退化为能量检测的变形价值会打折。4.2 核心函数Pietra-Ricci指数检测器下面这段Matlab函数是我一直在用的核心实现。function T pr_index_detector(X, noise_var) % X: Ns x Nt每行一个采样Nt个接收支路 % noise_var: 可选参数用于理论噪声特征值归一化 % T: PR指数检测统计量越大表示信号谱偏离越明显 [N, M] size(X); R (X * X) / N; % 样本协方差矩阵 lambda sort(eig(R), descend); % 特征值降序排列 % 归一化特征值除以均值消除噪声功率尺度影响 lambda_norm lambda / mean(lambda); % 理论参考CDF用单位均值指数分布近似噪声特征值分布(简化) % 若已知噪声功率可先把特征值除以noise_var再归一化 ref_cdf 1 - exp(-lambda_norm); % 经验CDF排序特征值的累积占比 obs_cdf (1:M) / M; % 计算Pietra-Ricci距离CDF差异的L1距离折半 T sum(abs(obs_cdf - ref_cdf)) / (2 * M); end这个实现用了单位均值指数分布作为理论参考严格来说不是精确的Marchenko-Pastur形式但好处是代码极短而且归一化后特征值确实围绕1波动PR指数的“相对变化”是稳健的。如果你想在论文里用严格版本可以把ref_cdf替换成按Marchenko-Pastur公式计算出的理论CDF。Marchenko-Pastur分布的参数 (c M/N)密度函数有闭式表达式Matlab写起来也不难我会在后面说明替换位置。4.3 多节点融合仿真主循环接下来是集中式数据融合的主仿真框架。这里采用“每个节点本地算PR指数融合中心做等增益合并”的结构。% 主仿真PR指数检测器 集中式数据融合 clear; clc; rng(2024); K 4; % 协作节点数 N 1000; % 每个节点采样点数 Nt 2; % 每个节点接收支路数 SNR_dB -20:2:0; MC 2000; Pfa_target 0.01; % 预分配 Pd zeros(numel(SNR_dB), 1); T_fusion_all zeros(MC, 1); for s 1:numel(SNR_dB) snr_lin 10^(SNR_dB(s)/10); for mc 1:MC % 生成主用户信号QPSK symbols (randi([0 1], N, 1) * 2 - 1) 1i * (randi([0 1], N, 1) * 2 - 1); sig_power var(symbols); noise_power sig_power / snr_lin; % 生成瑞利信道系数各节点独立 h (randn(K,1) 1i*randn(K,1)) / sqrt(2); % 各节点接收并计算本地PR指数 T_local zeros(K,1); for k 1:K noise sqrt(noise_power/2) * (randn(N,Nt) 1i*randn(N,Nt)); x h(k) * symbols * ones(1,Nt) noise; T_local(k) pr_index_detector(x, noise_power); end % 融合中心等增益合并 T_fusion mean(T_local); T_fusion_all(mc) T_fusion; end % 实际门限需要单独蒙特卡洛标定 % threshold quantile(T_noise_all, 1-Pfa_target); threshold 0.05; % 占位值正式仿真不可用 Pd(s) mean(T_fusion_all threshold); end figure; plot(SNR_dB, Pd, -o, LineWidth, 1.5); xlabel(SNR (dB)); ylabel(Detection Probability); grid on;注意上述主循环里threshold 0.05;是占位符实际跑的时候必须换成quantile(T_noise_all, 1-Pfa_target)。为了演示结构省略了标定段这在实际复现里要避免。一个干净的做法是写一个子函数calibrate_threshold(K,N,Nt,Pfa_target,MC)在纯噪声条件下统计融合后PR指数的分布并返回门限。4.4 结果怎么画看什么指标仿真结束后主要看两条指标一是检测概率 (P_d) 随信噪比变化曲线二是ROC曲线也就是 (P_d) 对 (P_{fa}) 的曲线。前者直观后者能体现检测器在全虚警区间的稳定性。把单节点K1和多节点K4或K8的曲线画在同一张图里就能看到协作增益。我自己跑出来的典型结果是目标虚警概率0.01采样点数1000二支路接收单节点在-8dB附近 (P_d) 到0.9四节点协作后同样 (P_d) 对应的SNR大约下移4~5dB八节点再下移1~2dB收益开始饱和。换句话说协作感知在前几个节点上收益最明显继续加节点边际递减。这一点在写报告或论文时值得重点分析不要只说“加了节点就变好”。5. 实战中踩过的坑与性能优化建议5.1 特征谱偏斜导致的误判PR指数对“特征谱整体形状”敏感但这也带来一个问题当采样点数特别少或接收矩阵近奇异时经验特征值分布会剧烈抖动。我遇到过在纯噪声情况下PR统计量也偏大的情况原因是仿真里噪声矩阵是复数高斯但忘了除2噪声功率被放大了一倍归一化又没用真实噪声功率做参考导致所有特征值都偏大理论CDF失配。排查办法是先把程序固定在纯噪声条件下打印归一化特征值的均值和方差。如果没有信号输入归一化特征值均值应该在1附近方差在可预测范围内如果均值严重偏离1说明归一化有问题。还有一种情况是数据矩阵里有NaN或者Inf特征值分解直接出错这类问题看起来像检测器不工作其实是数据预处理的问题。5.2 噪声功率未知场景的处理PR指数检测器不需要精确噪声功率但归一化方式不同会影响门限分布。我用的是“除以特征值均值”的归一化它天然不依赖噪声功率但代价是统计量和信号功率存在耦合。另一种做法是用一小段纯噪声参考样本估计噪声功率再对特征值归一化。前者简单后者更接近理论谱分布。如果场景里存在干扰信号或邻频泄漏特征谱可能被干扰峰污染PR指数容易把干扰当作主用户导致误警。这时可以在融合中心增加一个“特征值平滑”或者“去除最大特征值”的预处理再计算PR指数。不过这会降低对真实窄带信号的检测能力得不偿失需要根据目标主用户信号特点来定。5.3 融合权重不是越多越好集中式数据融合最容易犯的错误是“无脑等增益合并”。如果某个节点处于深度衰落或者上报信道质量很差它的本地统计量基本就是噪声直接参与均值合并会把整体统计量向H0方向拉低反而拖垮增益。更合理的做法是做加权合并权重可以是节点本地信噪比的估计值或者是节点历史检测可信度。哪怕只是最简单的置信度截断把低于某阈值的节点数据丢弃也比纯均值稳健。我实验过在四节点场景中把一个节点的信道故意设为极深度衰落等增益合并后性能比去掉该节点还差。后来改成信噪比加权合并每节点权重 (w_i \frac{1}{1\sigma_i^2/\hat{P}_i})性能才恢复到接近理想协作。这个细节在论文里不影响故事线但在工程项目里很关键。5.4 仿真耗时与代码加速技巧PR指数检测器要算特征值分解在蒙特卡洛大循环里会比较慢。我有几个加速技巧第一对固定维度的小矩阵直接用eig没问题但避免在循环里重复分配大数组第二可以把信号生成从“每节点循环”改成矩阵化生成一次性生成 (K * N) 的噪声矩阵再 reshape第三门限标定和性能测试分开标定用高精度门限比如20000次性能测试用5000次即可省时间第四如果只要检测趋势可以先跑一个低采样点数版本比如 N200跑通逻辑后再把 N 改到1000。Matlab的并行池也能加速蒙特卡洛循环但要注意随机数流一致性问题否则不同worker生成的噪声序列可能会导致统计结果波动。用parfor时最好在循环内显式设置独立的随机种子或者使用RandStream为每次迭代创建独立子流否则复现性很差。我用并行跑完一个多节点场景时间能压到原来的四分之一这是调参阶段的救命稻草。最后提一个容易被忽略的点如果你想把PR指数检测器和能量检测对比记得把两种检测器放到完全相同的门限标定流程下即都用蒙特卡洛标定各自的CFAR门限这样曲线差别才是真实算法差距而不是门限标定方法造成的假象。我以前吃过这个亏改完之后ROC曲线才真正拉开差距。
