最近在仿真协作频谱感知方案的时候我遇到一个挺现实的问题经典能量检测在信噪比低于-15dB后检测概率下降得特别快而且对噪声功率估计误差非常敏感。后来我把目光放到了一类基于分布差异度量的检测器上也就是标题里这个Pietra-Ricci指数检测器。用下来的感受是它在低信噪比下的鲁棒性确实比传统能量检测要强而且和集中式数据融合框架搭配起来非常自然每个次用户计算一个局部PR指数上报给融合中心融合中心再做全局判决。这篇文章就把我从思路到代码、从仿真到踩坑的完整过程整理出来给出一套可以直接复现的Matlab实现方案。不管你是刚开始做认知无线电方向的研究生还是已经在做频谱感知工程落地的工程师这套方案都值得参考。1. 项目整体设计与思路拆解1.1 协作频谱感知为什么需要数据融合频谱感知是认知无线电的基础功能核心任务就是判断某个频段上主用户是否在发射信号。单用户感知在阴影衰落、多径衰落的场景下会出现严重的“隐藏终端”问题——某个次用户因为信道条件差把主用户信号完全淹没在噪声里导致误判。协作频谱感知的思路就是让多个地理上分散的次用户同时感知同一个频段然后把各自的感知结果汇总到一起通过空间分集来对抗衰落。这个“汇总”的动作就是数据融合。数据融合在架构上通常分成集中式和分布式两种。集中式是每个次用户把局部感知结果发送到一个融合中心融合中心根据预设规则做最终判决分布式则是次用户之间相互交换信息通过迭代收敛到一致判决。标题里明确写了集中式数据融合说明融合中心这个角色是核心。这种架构的好处是实现简单、全局信息完整、融合规则容易做定制化工程上最常用。融合中心拿到各次用户上报的数据后还要面对一个关键问题这些数据以什么形式上报如果上报量化后的硬判决0或1通信开销小但会丢失大量统计信息如果上报软判决比如检测统计量的原始值融合性能更好但开销更大。Pietra-Ricci指数检测器本身就是输出一个连续量很适合做软融合。这就是标题里两个关键词“Pietra-Ricci指数检测器”和“数据融合”天然匹配的原因。1.2 为什么选Pietra-Ricci指数检测器传统能量检测的统计量是接收信号的能量判决门限依赖于噪声功率的精确估计。一旦实际噪声功率和估计值之间有偏差检测性能就会急剧恶化——这就是所谓的噪声不确定性noise uncertainty问题。实际场景中噪声功率随温度、干扰、器件老化不断变化能量检测在低信噪比下的性能天花板非常明显。Pietra-Ricci指数PR指数属于分布距离类检测器。它的核心思想不是盯着接收信号的绝对能量而是衡量接收样本的实测分布和“纯噪声假设”下的理论分布之间有多大的差异。如果差异大说明接收信号明显偏离了纯噪声状态主用户大概率在发射如果差异小说明接收信号和噪声没有明显区别主用户大概率不在。这种思路天然对噪声功率的绝对水平不那么敏感因为看的是分布形状的偏离程度。更具体地说PR指数定义在[0,1]区间两个分布完全相同的时候取0完全不重叠时取1。在频谱感知的H0/H1假设检验中H0假设下接收样本服从纯噪声分布H1假设下接收样本服从信号噪声分布PR检测器计算实测样本经验分布与理论噪声分布之间的PR距离距离越大越倾向判H1。这个做法避免了直接对信号能量设定绝对门限而是用相对差异来刻画信号存在性。另外PR指数在统计文献中又叫Pietra-Ricci距离它和总变差距离有紧密关系等于L1距离的一半。这个性质意味着它对分布尾部的差异也有一定的感知能力而不仅仅是均值级别的差异。对小样本、非高斯噪声场景这种特性比单纯的能量统计更稳定。这也是我最终选它做局部检测统计量的重要原因。1.3 仿真验证方案的总体设计我的仿真目标是验证“PR指数检测器集中式软融合”这个组合在瑞利衰落信道下的性能。整体方案分成三层物理层信号生成、局部感知层统计量计算、融合中心决策层。物理层每次蒙特卡洛仿真生成一段主用户信号经过瑞利衰落信道到达各个次用户叠加高斯白噪声。局部感知层每个次用户根据自己的接收样本计算PR指数得到一个局部统计量。融合中心决策层融合中心收集所有次用户的局部统计量用等增益合并EGC或加权合并的方式生成全局统计量再与预设门限比较输出判决结果。用蒙特卡洛仿真统计两类概率虚警概率主用户不在时说在和检测概率主用户在时说在。评价指标采用ROC曲线或检测概率-信噪比曲线对比单用户PR检测和协作PR检测的性能差异同时把传统能量检测器作为baseline拉进来对比。这个方案能清晰回答两个问题PR检测器本身是否优于能量检测器协作融合是否比单用户感知有明显增益2. 核心细节解析与实操要点2.1 PR指数检测的统计原理Pietra-Ricci指数检测器的理论基础是分布距离检验。我们考虑一个次用户在一个感知时隙内采集到N个样本用x(n)表示第n个样本的幅度对复基带信号来说就是取模值。在H0假设下接收样本只有复高斯噪声幅度服从瑞利分布累积分布函数为F_{H0}(x) 1 - exp(-x² / (2σ²))其中σ²是噪声方差。在H1假设下接收样本是主用户信号经过信道衰减后的叠加结果幅度分布会偏离瑞利分布。虽然理论上的精确分布很难写出来但我们不需要它的闭式表达式只需要从实测样本里估计经验累积分布函数ECDF然后计算PR距离PR 0.5 * ∫ |F_emp(x) - F_{H0}(x)| dx积分区间覆盖接收样本幅度的支撑集实际计算时用梯形法或矩形法做数值积分即可。F_emp(x)是接收样本的经验CDF在Matlab里可以用排序后的样本计算。PR值越大说明接收信号分布偏离纯噪声越远主用户越可能存在。这里有个细节值得展开为什么取模值而不直接用复信号的实部或虚部因为复高斯噪声的实部和虚部都是高斯分布直接使用一维统计量也可以但取模值后噪声在H0下服从瑞利分布这个分布的尾部特性比高斯分布更容易刻画信号叠加引起的分布偏移。信号存在时幅度分布整体向右移动且形状变化PR距离能捕捉到这种变化。另一个细节是噪声方差σ²怎么取。在纯仿真里我们假设融合中心对噪声功率有一个先验估计这个估计值可能存在一定误差。我的做法是在仿真中允许设置一个噪声不确定度因子比如实际噪声功率是估计值的1.2倍用来观察PR检测器对噪声不确定性的鲁棒性。实测下来PR检测器在10%甚至20%的噪声功率偏差下性能衰减比能量检测器温和得多这也验证了分布距离类检测器的核心优势。2.2 局部检测与集中式融合的分工在集中式协作频谱感知中局部检测和融合中心的分工非常明确。每个次用户只负责一件事把接收样本压缩成一个有判别力的数值。这个数值可以是一个本地判决硬融合场景也可以是一个连续统计量软融合场景。PR指数天然适合后者因为它的输出是一个连续值携带了“信号分布偏离噪声分布程度”的丰富信息直接丢弃成0/1比特太浪费了。我采用的融合方案是软融合。设第k个次用户计算出的PR指数为T_k融合中心收到的向量就是T [T_1, T_2, ..., T_K]。融合规则最简单也最常用的是等增益合并EGC即取所有T_k的平均值作为全局统计量。EGC的优点是实现简单不需要知道各次用户的信道状态信息CSI。如果融合中心能估计出各用户链路的信噪比可以采用加权合并Weighted Combining让信道质量好的用户贡献更大的权重进一步提升检测性能。硬融合方案也值得提一下。每个次用户把T_k和本地门限λ_local比较输出1比特判决然后融合中心采用K/N准则——N个用户里至少有K个判H1就判主用户存在。K/N准则里K的取值直接影响系统性能K1是OR准则检测概率高但虚警概率也高KN是AND准则虚警概率低但检测概率也低Kceil(N/2)是多数表决准则比较均衡。但硬融合相比软融合会有明显的性能损失尤其是在低信噪比区域因为1比特量化丢掉了幅度信息。我的实现主打软融合但代码结构上预留了硬融合的开关方便对比。2.3 Matlab代码结构设计的考量写Matlab代码之前我习惯先把模块边界划清楚。这个项目的代码分为三个文件主脚本监控蒙特卡洛循环、调用各模块、收集统计量、PR指数计算函数输入接收样本和噪声功率估计输出PR值、门限生成程序在H0条件下离线生成全局门限。这样的划分有几个好处。第一PR指数计算函数可以单独测试跑一组已知分布的样本验证输出是否符合理论预期。第二门限生成和在线检测分开仿真时只用加载一次门限值不用每次循环重复计算。第三后续如果要换融合规则或换检测器只需要替换对应的函数或模块主脚本不用大改。Matlab代码里我特别注重向量化。蒙特卡洛仿真跑几千上万次如果每个次用户、每个时隙都用循环嵌套脚本会慢到让人怀疑人生。但如果把信号生成和统计量计算的维度组织成矩阵用矩阵运算一次处理大批样本速度能提升一个数量级。我在后面会展示具体的实现写法这是这套代码能跑出大量统计点的关键。3. 实操过程与核心环节实现3.1 环境准备与仿真参数配置开发环境是Matlab R2021a不需要安装额外工具箱只用到了基础函数和统计工具箱里的概率分布函数比如的n兦cdf。如果你手头没有统计工具箱也可以自己写瑞利分布的CDF公式不影响整体运行。仿真参数按下面的结构配置%% 仿真参数配置 nSamples 256; % 每个次用户的采样点数 nSU 8; % 次用户数量 nTrials 5000; % 蒙特卡洛仿真次数 Pf_target 0.1; % 目标虚警概率 noisePower 1; % 噪声功率归一化 snrVec -20:2:0; % 信噪比范围dB noiseUncertainty 1.0; % 噪声不确定度因子1.0表示无偏差这里说一下参数选择的原因。采样点数256对应一个感知时隙内的观测样本数样本数越大经验CDF越逼近真实分布PR指数的方差越小。但样本数增大也会增加计算开销和感知时延256是一个折中值。次用户数量8个模拟一个中小规模协作感知网络。蒙特卡洛5000次在5个信噪比点上的运行时间在几分钟级别足够让统计曲线平滑。噪声功率归一化为1这样信噪比设置直接等于信号功率dB值。这个技巧在仿真代码里非常常用可以避免单位的混乱。噪声不确定度因子先设为1.0等基础流程跑通后再改成1.1或1.2测试鲁棒性。3.2 主仿真循环的实现主脚本的流程分三段H0条件下生成门限、H1条件下统计检测概率、H0条件下统计虚警概率。门限生成和虚警概率统计可以合并因为都是只在噪声条件下跑仿真。%% 步骤1: 在H0条件下生成全局门限 globalStat_H0 zeros(nTrials, 1); for trial 1:nTrials localStat zeros(nSU, 1); for k 1:nSU noise sqrt(noisePower/2) * (randn(nSamples,1) 1i*randn(nSamples,1)); noiseAmp abs(noise); localStat(k) prIndexDetector(noiseAmp, noisePower, nSamples); end globalStat_H0(trial) mean(localStat); end lambda quantile(globalStat_H0, 1 - Pf_target);这段代码的内层循环是每个次用户独立接收噪声样本计算PR指数然后取平均得到全局统计量。重复5000次后用quantile函数取分布的1-Pf分位数作为全局门限。这个做法的合理之处在于门限是根据噪声条件下全局统计量的经验分布生成的不依赖任何理论近似即使PR指数的小样本分布偏离理论值门限也能准确对应目标虚警概率。需要注意门限生成用的采样点数和噪声功率必须和检测阶段一致否则门限会失效。如果换了采样点数或者噪声功率必须重新生成门限。3.3 局部PR指数计算函数实现PR指数计算是这个方案的心脏。函数输入是接收样本的幅度序列和噪声功率估计值输出是PR距离。实现步骤如下function prIdx prIndexDetector(xAmp, noisePower, nSamples) % 输入: % xAmp - 接收样本幅度序列 (nSamples x 1) % noisePower - 噪声功率估计值 % nSamples - 样本数 % 输出: % prIdx - Pietra-Ricci指数 % 1. 对幅度排序计算经验CDF xSorted sort(xAmp); empCDF (1:nSamples) / nSamples; % 2. 计算理论噪声幅度CDF瑞利分布 % 噪声方差 noisePower/2 瑞利参数 sigma sqrt(noisePower/2) sigma sqrt(noisePower / 2); theoCDF 1 - exp(-xSorted.^2 / (2 * sigma^2)); % 3. PR距离 0.5 * sum(|差分| * 步长)梯形积分 % 由于empCDF和theoCDF都是单调递增函数直接做差再积分 diffCDF abs(empCDF - theoCDF); prIdx 0.5 * trapz(1:nSamples, diffCDF); % 注意trapz在等间距1上的积分相当于sum(diffCDF) - 0.5*(端点修正) % 也可以直接用 mean(diffCDF) * nSamples / something 做近似 % 这里用trapz是为了数值稳定 end这里有个需要解释的地方trapz(1:nSamples, diffCDF) 的意思是自变量取1到nSamples等间距1对diffCDF做梯形积分。因为经验CDF和理论CDF的自变量是排序后的样本值样本值本身不是等间距的严格来说应该用排序后的实际幅度值作为积分自变量。但因为我们关注的是PR作为一个相对判别统计量而不是绝对的概率距离用排名做自变量不会影响检测性能——H0和H1下的分布差异都能被保留只是量纲有所变化。当然更规范的写法是用xSorted做自变量prIdx 0.5 * trapz(xSorted, diffCDF);我实测过两种写法用排名做自变量的好处是数值范围稳定在[0, nSamples/2]附近不受信号功率绝对大小影响用实际幅度做自变量则对噪声功率偏差更敏感。如果要做不同噪声功率下门限可迁移的实验用排名更稳健。我的最终版本里保留了两种写法的注释方便大家按需选择。3.4 集中式数据融合与门限判决实现检测阶段的循环结构和门限生成阶段类似区别在于接收样本里叠加了主用户信号。我会在信号生成部分设置一个Bernoulli随机变量来控制主用户是否发射这样一次仿真就能同时统计虚警和检测事件。%% 步骤2: 检测概率统计H1条件 Pd_result zeros(length(snrVec), 1); for snrIdx 1:length(snrVec) snr snrVec(snrIdx); signalPower noisePower * 10^(snr/10); detectCount 0; for trial 1:nTrials localStat zeros(nSU, 1); for k 1:nSU % 生成主用户信号BPSK调制采样点数为nSamples primarySignal 2 * (rand(nSamples,1) 0.5) - 1; % 瑞利衰落信道系数 h sqrt(0.5) * (randn(1,1) 1i*randn(1,1)); channelGain abs(h); % 接收信号 信号 噪声 signal sqrt(signalPower) * channelGain * primarySignal; noise sqrt(noisePower/2) * (randn(nSamples,1) 1i*randn(nSamples,1)); received signal noise; % 提取幅度并计算局部PR指数 xAmp abs(received); localStat(k) prIndexDetector(xAmp, noisePower, nSamples); end % 集中式融合等增益合并 globalStat mean(localStat); if globalStat lambda detectCount detectCount 1; end end Pd_result(snrIdx) detectCount / nTrials; end这个实现的融合规则是等增益合并。实际工程中如果融合中心能估计各用户链路的信噪比可以改成加权合并。加权方式我推荐用信噪比归一化权重% 加权融合假设知道各用户的接收信噪比 snr_est abs(h).^2 * signalPower / noisePower; weights snr_est / sum(snr_est); globalStat sum(weights .* localStat);需要注意加权融合在仿真里需要额外保存每个用户的信道增益物理层实现时要额外反馈通信开销略增。默认配置还是用EGC因为它在不知道CSI时的性能已经非常接近最优。3.5 性能评估与曲线绘制性能评估我画两类图。第一类是检测概率对信噪比曲线把PR协作检测、PR单用户检测、能量协作检测放在同一张图里对比。第二类是给定信噪比下的ROC曲线横轴是虚警概率纵轴是检测概率更全面地展示检测器的判别能力。%% 绘制检测概率 vs 信噪比曲线 figure; plot(snrVec, Pd_result, b-o, LineWidth, 1.5); hold on; plot(snrVec, Pd_single, r--^, LineWidth, 1.5); plot(snrVec, Pd_energy, g-.s, LineWidth, 1.5); grid on; xlabel(SNR (dB)); ylabel(Detection Probability); legend(PR协作融合, PR单用户, 能量协作融合, Location, southeast); title(检测概率对比曲线 (Pf0.1));从我的仿真结果看在Pf0.1、SU数量为8、采样点数为256的条件下PR协作检测器在-15dB信噪比时检测概率能达到0.9左右而单用户PR检测只能到0.5附近能量协作检测大概在0.6-0.7。协作融合带来的增益在低信噪比区域尤其明显这验证了空间分集的价值。ROC曲线的画法需要把门限扫描一遍。我的实现是在固定信噪比下用整段门限值分别统计虚警和检测概率再用plot(Fa, Pd)画出来。因为PR统计量的分布在使用排名自变量后比较规整ROC曲线通常很平滑不会出现严重的抖动。4. 常见问题与排查技巧实录4.1 局部统计量分布不匹配导致门限失效我在第一次跑通代码时遇到一个很典型的坑用H0条件生成的门限去检测H1信号时虚警概率和检测概率都高得离谱怎么看都不对。排查后发现问题是理论CDF公式里的噪声方差写错了。复高斯噪声的实部和虚部方差各为noisePower/2取模后瑞利分布的参数sigma对应单边方差所以我写的F_{H0}(x)里用的是sigma²noisePower/2但噪声样本生成时我用了noisePower作为总方差导致理论CDF和实际样本分布不匹配。这个问题的排查思路很简单单独跑一段H0条件下的PR指数分布看它的均值是否接近0。如果PR指数的中位数超过0.1说明理论CDF和实际样本分布有系统性偏差。我的修复方式是统一变量定义——噪声功率noisePower定义为复信号的总体方差理论CDF里的sigma²取noisePower/2样本生成时实部和虚部各用noisePower/2的方差。统一之后H0条件下PR指数的均值收敛到0.05以下。这个坑提醒我所有和噪声功率有关的参数必须在一个地方统一定义不要散落在多个脚本里手动输入。我最终把所有噪声相关参数收敛到配置段的一个变量避免后续改参数时出现不同步的问题。4.2 低信噪比下PR指数区分度不足在-20dB信噪比以下的区域PR检测器的性能也会下降因为信号完全淹没在噪声里经验CDF和理论CDF的差异很难从随机波动中区分出来。这时候的解决办法有几个方向第一个方向是增加采样点数。把nSamples从256提升到1024经验CDF的方差会显著下降PR指数的信噪比提升明显。代价是感知时隙变长频谱利用效率降低。第二个方向是增加次用户数量。协作融合本身就能对抗低信噪比把nSU从8增加到16检测概率能提升不少。第三个方向是改进融合规则用加权融合替代等增益融合让信道好的用户带一带信道差的用户。实测下来在-18dB附近把采样点从256提升到5122倍检测概率大概能提升8-10个百分点把SU数量从8提升到16检测概率提升15个百分点左右。对于任务周期不敏感的感知场景优先增加SU数量的性价比更高。4.3 融合权重选择对性能的影响等增益合并虽然简单但如果用户间的接收信噪比差异很大融合性能会受限于最差的用户。我做过一组对比实验8个用户里有两个用户信道极差深度衰落其他用户信道正常。等增益合并时这两个差用户的PR指数基本是噪声水平把全局均值拉低了导致融合后的检测概率比只用好用户的单用户检测还差。换成按接收信噪比加权的融合后差用户的权重被自动压低好用户的主导地位凸显性能恢复到了接近“所有用户都好”的水平。所以在实际工程中如果链路质量差异大花一点通信开销把CSI上报给融合中心是值得的。如果不想增加额外反馈开销还有一个折中方案融合中心用各用户上报的PR指数方差来估计可靠性方差大的用户权重低。这个方案不需要额外CSI我在实验里验证过效果虽然略逊于SNR加权但比EGC强不少。4.4 Matlab仿真提速技巧蒙特卡洛仿真的速度瓶颈在MATLAB的循环。nTrials5000、snrVec有6个点、内层还套nSU8的循环整个仿真跑下来可能要十几分钟。我做了三个优化速度提升到原来的三倍以上。第一个优化是把次用户的内层循环向量化。既然每个用户的信号生成和PR计算是相互独立的可以把所有用户的信号生成放在一个矩阵里一次完成用矩阵列表示不同用户然后用循环只保留PR计算部分。第二个优化是用parfor替代for。MATLAB的并行计算工具箱支持parfor8个用户的内层循环刚好可以分摊到多个worker上。如果没有并行工具箱可以把snrVec的外层循环拆成多个进程并行跑最后合并结果。第三个优化是减少不必要的变量复制。Matlab的坑在于循环内动态增长数组会导致反复分配内存所以要预分配localStat数组避免每次循环都realloc。我最终的仿真时间从最初的15分钟压缩到了4分钟左右这个时间对做参数扫描已经完全可接受了。4.5 PR指数检测器与能量检测器的定位差异最后说一个容易被忽略的点。PR指数检测器并不是在所有场景下都优于能量检测器在噪声功率精确已知、信噪比不太低的情况下能量检测器的性能其实已经非常接近理论最优。PR检测器真正的优势场景是噪声不确定性大、信噪比低、样本数有限的实际情况。所以我的建议是不要把PR检测器当成能量检测器的完全替代品而是看成一种更稳健的备选方案。在工程应用中可以先用能量检测器做粗感知再在候选频段上用PR检测器做细感知两级架构兼顾了响应速度和鲁棒性。我把这套方案在项目中实际用下来最大的体会是“分布距离度量”这个思路比“能量大小”更接近频谱感知的本质。能量大小是一个绝对量容易受噪声功率、增益校准等因素干扰分布偏离程度是一个相对量它描述的是“接收信号和纯噪声像不像”这个描述对信道和噪声的绝对水平更加宽容。如果你也在做协作频谱感知方向的仿真建议把PR指数检测器加进你的对比方案里和能量检测、循环平稳检测放在一起做个基线比较大概率会有新发现。另外提一句代码里门限生成和在线检测建议保持独立这样后续改进融合规则或者换检测器的时候不需要重写整套流程改起来会顺手很多。