基于Pietra-Ricci指数检测器的协作频谱感知及Matlab实现-嵌云网-嵌入式AI开发资源站

基于Pietra-Ricci指数检测器的协作频谱感知及Matlab实现

小糖元·2026-09-13·阅读 21 GEO检测（免费）

协作频谱感知Pietra-Ricci指数检测器集中式数据融合

|难度：普通|阅读约 23 分钟

做认知无线电的朋友对协作频谱感知应该都不陌生。多节点把感知结果汇总到融合中心，确实比单节点硬扛衰落和阴影可靠得多。不过网上大多数教程翻来覆去都是能量检测加硬判决融合，看多了总觉得缺了点新东西。这次我分享一个不太一样的方案：用Pietra-Ricci指数检测器，在集中式数据融合架构下做协作频谱感知，并给出完整的Matlab实现思路。

这个方法的核心不是比谁的信号能量高，而是比较“接收样本的经验分布”和“纯噪声的理论分布”差多少。只要有主用户的信号混进来，样本分布就会偏离噪声分布，这个偏离程度用Pietra-Ricci指数量化，融合中心拿这个指数做判决。对于想尝试非参数检测、或者做协作频谱感知算法对比的朋友，这个思路值得收藏。

先搞清楚Pietra-Ricci指数检测器是什么

1.1 把频谱感知当成“分布差异检测”

频谱感知在数学上就是一个二元假设检验问题：

H0：接收信号里只有噪声，r(n) = w(n)
H1：接收信号里包含主用户信号，r(n) = s(n) + w(n)

传统能量检测的做法是计算接收信号的能量，然后跟一个能量门限比较。这个方法简单，但门限高度依赖噪声功率的估计，一旦噪声功率估不准，检测性能会急剧恶化。

Pietra-Ricci指数检测器换了一个思路：不直接算能量，而是先构建接收样本的经验累积分布函数，再跟纯噪声的理论累积分布函数做比较。如果实际环境是H0，经验分布应该和理论噪声分布很接近；如果是H1，样本里混入了主用户信号，分布就会发生形变，两个CDF之间的差异会变大。

工程上我采用的Pietra-Ricci统计量定义如下：

D = ∫ |F_N(t) - F_0(t)| dt

其中F_N(t)是接收样本的经验CDF，F_0(t)是噪声的理论CDF。D越大，说明接收样本越不像是纯噪声，越倾向判为“主用户存在”。

你可能想问，这个和KL散度、最大均值差异（MMD）那类分布距离有什么区别？区别在于Pietra-Ricci指数是基于CDF的L1距离，计算简单，对尾部分布的敏感度适中，不需要估计概率密度函数，直接用排序样本做梯形积分就能得到，非常适合Matlab实现。

1.2 集中式数据融合在这个方案中的定位

协作频谱感知有两种典型架构：

分布式融合：节点之间相互交换本地判决或中间统计量，各自形成最终判决；
集中式融合：所有节点把本地观测数据或统计量上传到融合中心，由融合中心统一处理并给出全局判决。

本文用的是集中式数据融合。各次用户（SU）在感知时隙内采集N个复基带采样，然后把原始采样数据（或经过功率归一化的数据）送到融合中心（FC）。FC把所有节点的样本拼接成一个大的样本集合，再计算Pietra-Ricci统计量。

为什么集中式比分布式更适合这个检测器？因为Pietra-Ricci指数需要足够多的样本才能把经验CDF做得平滑。单个节点N个样本，分布估计误差大；K个节点都上传后，融合中心手里就有K×N个样本，CDF逼近真实分布的能力强很多，检测统计量也更稳定。简单说，集中式融合是用通信开销换检测性能。

系统模型与核心参数设计

2.1 网络模型与信道假设

先搭一个仿真模型。假设系统里有：

一个主用户（PU），发射复基带信号s(n)，调制方式不限，常见用QPSK或OFDM；
K个次用户（SU），每个SU在感知时隙内采集N个采样；
一个融合中心（FC），负责收集K个节点的采样并做集中式判决；
每个SU到PU之间的信道是独立瑞利衰落，再加复高斯白噪声。

第k个节点接收信号可以写为：

r_k(n) = h_k × s(n) + w_k(n)

其中h_k是复信道增益，w_k(n)是循环对称复高斯噪声，实部和虚部方差各为σ²。为了简化，仿真时可以假设报告信道（节点到FC）是无差错、无延时的理想信道，先专注于感知算法本身。

2.2 检验统计量到底怎么构造

集中式融合中心拿到K个节点、每个N个采样后，把所有样本排列成一个长度为M = K × N的向量r_all。接下来有两种常见处理方式：

直接使用复采样值：经验CDF和理论CDF都基于复数的实部或虚部，但需要分别处理；
使用能量样本：对每个采样取模平方 y = |r|²，在H0下y服从指数分布，理论CDF有闭式表达式。

我在工程实现中更喜欢第二种。原因是模平方同时利用了实部和虚部的信息，检测信息不浪费，而且理论CDF非常简单。

假设噪声复方差为σ_w²，也就是实部虚部方差各为σ_w²/2，那么H0下y = |w|²的CDF是：

F_0(t) = 1 - exp(-t / σ_w²)

注意这里分母是σ_w²，不是2σ_w²。原因在于复高斯噪声w的实部和虚部方差都是σ_w²/2，所以|w|²的均值等于σ_w²。

经验CDF定义为：

F_M(t) = (1/M) × Σ I(y_i ≤ t)

然后统计量：

D = ∫₀^∞ |F_M(t) - F_0(t)| dt

D的物理意义很直观：它衡量了观测能量样本分布和纯噪声能量分布之间的总偏差面积。

2.3 门限设定与蒙特卡洛仿真思路

D的解析分布很难推导，所以工程上基本都用蒙特卡洛仿真来定门限。具体做法是：

在H0条件下生成大量噪声样本，计算对应的D0；
对D0做经验排序，取第(1 - P_fa)分位点作为检测门限γ；
在H1条件下生成含信号样本，计算D1；
若D1 > γ，判为H1，统计检测概率P_d。

需要注意，门限只和噪声分布、样本量有关，和信噪比无关，所以同一组门限可以用于所有SNR点的检测概率计算，不用每组SNR重新算门限。这个细节能省不少仿真时间。

Matlab代码实现与模块拆解

3.1 主程序流程

先给一个主程序框架，方便你直接照着搭。

matlab复制clearvars; close all; clc;
rng(2024);

% 系统参数
K = 4;                  % 次用户数
N = 1024;               % 每个节点采样点数
numMC = 2000;           % 蒙特卡洛次数
Pfa_target = 0.1;       % 目标虚警概率
SNR_dB = -14:2:0;       % 信噪比范围

sigma_w2 = 1;           % 复噪声功率，实虚部方差各为0.5

% 预分配
D0 = zeros(numMC, 1);
Pd = zeros(length(SNR_dB), 1);

% 第一步：H0门限估计
fprintf('Estimating threshold under H0...\n');
for mc = 1:numMC
    r = sqrt(sigma_w2/2) * (randn(K, N) + 1j*randn(K, N));
    D0(mc) = calcPietraRicci(r(:), sigma_w2);
end
threshold = quantile(D0, 1 - Pfa_target);

% 第二步：H1检测概率
for snrIdx = 1:length(SNR_dB)
    SNR_lin = 10^(SNR_dB(snrIdx)/10);
    Ps = SNR_lin * sigma_w2;   % 信号功率
    detCount = 0;

    for mc = 1:numMC
        % 主用户信号，QPSK
        s = (2*randi([0 1], 1, N) - 1 + 1j*(2*randi([0 1], 1, N) - 1)) / sqrt(2) * sqrt(Ps);
        
        % 多节点接收
        r_all = zeros(K, N);
        for k = 1:K
            h = (randn + 1j*randn) / sqrt(2);   % 瑞利衰落
            r_all(k, :) = h .* s + sqrt(sigma_w2/2) * (randn(1, N) + 1j*randn(1, N));
        end
        
        % 融合中心计算统计量
        D = calcPietraRicci(r_all(:), sigma_w2);
        if D > threshold
            detCount = detCount + 1;
        end
    end
    Pd(snrIdx) = detCount / numMC;
    fprintf('SNR = %.1f dB, Pf = %.2f, Pd = %.4f\n', SNR_dB(snrIdx), Pfa_target, Pd(snrIdx));
end

3.2 核心函数：计算Pietra-Ricci统计量

这是整个实现的重头戏。需要特别注意积分边界问题。

matlab复制function D = calcPietraRicci(r, sigma_w2)
    % r: 所有节点的复采样样本，列向量；sigma_w2: 复噪声功率
    y = abs(r(:)).^2;               % 能量样本
    y = sort(y);
    M = length(y);

    % 经验CDF在样本点的值
    Fn = (1:M)' / M;

    % 理论噪声CDF：1 - exp(-t / sigma_w2)
    F0 = 1 - exp(-y / sigma_w2);

    % 主体积分：用梯形法
    D = trapz(y, abs(Fn - F0));

    % 第一段：t从0到y(1)，此时Fn=0
    D = D + (y(1) - sigma_w2 * (1 - exp(-y(1)/sigma_w2)));

    % 尾段：t从y(M)到无穷，此时Fn=1，F0趋向于1
    D = D + sigma_w2 * exp(-y(M)/sigma_w2);
end

为什么这段代码能提升精度？如果只做trapz，实际上忽略了0到第一个样本点之间、以及最后一个样本点往后的分布差异。当样本量M很大时这两个边界项可以忽略，但M只有几百的时候，边界项对D的贡献很大，直接影响检测门限，所以必须补上。

3.3 融合中心的样本拼接与判决

融合中心的功能就是把所有节点样本拼成一个长向量，再调用上面的核心函数。我在主程序里用r_all(:)这一步完成了拼接，效果等同于：

matlab复制r_combined = [];
for k = 1:K
    r_combined = [r_combined; r_all(k, :)'];
end
D = calcPietraRicci(r_combined, sigma_w2);

要注意样本拼接顺序不影响经验CDF，因为排序后顺序无关。融合中心不需要知道哪个样本来自哪个节点，只需要统一做分布比较。这也是集中式数据融合的一个优势。

3.4 结果可视化和能量检测对比

画一条检测概率曲线：

matlab复制figure;
semilogy(SNR_dB, max(Pd, 1e-3), 'bo-', 'LineWidth', 1.5);
grid on;
xlabel('SNR (dB)');
ylabel('Detection Probability');
title('Pietra-Ricci Detector for Cooperative Sensing');

如果你想和能量检测对比，可以在同样的系统模型下实现一个集中式能量检测器：融合中心把所有能量样本平均，得到平均能量T，门限也通过H0蒙特卡洛仿真得到，最后同样统计Pd。两个曲线画在同一张图里，就能直观看出差异。

实验结果与关键参数影响

4.1 不同SNR下的检测性能

以K=4、N=1024为例，在我这组参数下测试，检测概率随SNR的变化呈现出典型的“S型”曲线。低信噪比区域，比如-14dB到-10dB，检测概率上升缓慢；到了-6dB以后，曲线快速拉起；接近0dB时，检测概率基本逼近1。

表格形式做个示意（具体数值随噪声和信道随机种子变化）：

SNR (dB)	检测概率Pd
-14	0.07
-12	0.13
-10	0.24
-8	0.46
-6	0.73
-4	0.92
-2	0.99

这个曲线的特点是：低SNR区域没有能量检测那么激进，一旦越过某个阈值，检测概率增长非常快。在协作场景下，节点数越多，曲线斜率越陡，低SNR性能提升越明显。

4.2 节点数K和采样点数N怎么影响结果

集中式Pietra-Ricci检测器最吃的是样本量。M = K × N直接决定经验CDF的平滑程度。如果M太小，经验CDF毛糙，统计量D的波动大，门限和检测概率都不可靠。

从实际仿真来看，K=2、N=512时，低SNR下检测性能一般；K=4、N=1024时，曲线明显改善；继续加到K=8、N=2048，曲线的提升幅度会变小，但仿真时间成倍增加。所以建议先从K=4、N=1024起步，跑通了再调参数。

节点数K和采样点数N对系统的影响不完全相同：

K增大，相当于融合了更多独立信道的观测，对抗衰落和阴影的能力更强；
N增大，只是让单个节点对噪声分布的估计更准，对分布式协作增益没有直接贡献。

所以协作频谱感知里通常先加节点，再考虑增加每节点采样长度。

4.3 与能量检测器的性能对比

能量检测器在高SNR下表现不差，但在噪声不确定性存在时非常脆弱。Pietra-Ricci指数检测器的一个明显优势是它对信号调制方式不敏感。只要主用户信号改变了接收样本的能量分布，经验CDF和噪声CDF之间就会出现面积偏差，检测器就能捕获到。

另一个优势是它对非高斯噪声的鲁棒性更好。能量检测只关注二阶矩，而非高斯噪声的二阶矩往往不稳定，门限容易失真。Pietra-Ricci检测器比较的是整个分布形状，对某些重尾噪声的适应能力更强。

当然它也有代价。计算量比能量检测大，需要排序和梯形积分；另外它需要相对准确的噪声分布模型，如果噪声方差估计偏差较大，理论CDF本身就不准，检测器性能同样会下降。

常见问题与排查技巧实录

5.1 统计量数值计算总是偏小

这个坑我一开始踩过。最开始只用trapz计算主体积分，忽略了首段和尾段，结果门限估出来明显偏低，H0下虚警率远超设定值。后来把边界项补上才恢复正常。

另一个可能导致偏小的原因是y排序后跨度太大，trapz的网格不均匀。解决办法是不用均匀网格，直接用排序后的样本点作为积分节点。这样经验CDF的突变点都在积分节点上，梯形积分误差可控。

5.2 门限仿真结果不稳定

蒙特卡洛次数太少时，门限会抖得很厉害。建议numMC至少2000，如果计算资源允许，最好到5000以上。门限估计和检测概率仿真要用独立的噪声样本集，不然会产生人为相关性，导致曲线偏乐观。

还有一个容易忽略的点：复噪声功率sigma_w2的单位。很多同学把复高斯噪声写成randn + 1j*randn，然后默认噪声功率为1，但实际上这个时候实部虚部各有一个单位方差，总功率是2。我在代码里用sqrt(sigma_w2/2)生成实部和虚部，保证总噪声功率等于sigma_w2。如果这里写错，理论CDF里的指数分布参数也会跟着错，整个检测器就废了。

5.3 仿真速度太慢怎么办

三层循环（SNR、蒙特卡洛、节点）跑起来确实慢。可以先做两件优化：

把信号生成、信道衰落和噪声生成尽量向量化，避免内层节点循环；
如果安装了Parallel Computing Toolbox，把蒙特卡洛循环改成parfor。

我实测parfor在四核机器上能快2到3倍。注意parfor里的随机数生成要用不同的随机流，避免每个worker产生相同序列。最简单的方法是在循环开始时用RandStream设置不同子流，或者接受默认全局流的分段行为。

5.4 画ROC曲线时出现“反曲线”

如果检测概率在低虚警概率时反而低于高虚警概率，通常是门限扫描方向写反了。Pietra-Ricci检测器判决规则是“D > γ判H1”，所以门限γ从小到大扫，虚警概率从高到低。画图时横坐标P_fa应该用1 - P_fa的分位点去对应。

还有一个情况：如果噪声方差设置不当，H0和H1的分布差异可能反号，这时统计量不仅不增大，反而可能减小。遇到这种情况，先单独打印H0和H1下的D分布，确认H1均值和分位数确实高于H0，再做后续分析。

一些项目落地时的个人体会

这套代码跑下来，我最深的感受是：Pietra-Ricci检测器不是一个“无脑强”的算法，但它给了频谱感知一个非常灵活的建模角度。你不需要去猜主用户信号的调制格式，也不需要精确估计信号功率，只需要对噪声分布建模，然后把“异常”交给分布差异去捕捉。这种思路在信号特征不明、干扰来源复杂的场景下尤其有价值。

在实际使用中，我一般会先跑一遍H0下的D分布，看看统计量的数值范围，再决定门限和样本量。这个步骤看起来很基础，却能避免后面很多返工。另外，如果你想把算法推广到宽带频谱感知，可以考虑对每个子带分别构造Pietra-Ricci统计量，再做融合，思路是完全相通的。

最后分享一个小技巧：在写Matlab代码时，尽量把检测统计量封装成独立函数，不要把所有逻辑都堆在主脚本里。这样后续换数据、换噪声模型、换融合策略都只需要改一小部分代码。我这套代码里的calcPietraRicci函数就是改了三版才稳定下来，前期一直把积分逻辑嵌在主循环里，调试一次要跑半天，封装之后效率高多了。