1. 协作频谱感知与集中式数据融合为什么要抱团取暖说句实话做频谱感知最郁闷的不是算法写不出来而是单节点在低信噪比环境里怎么调都压不住虚警。我最早做认知无线电仿真的时候一股脑用能量检测结果在阴影衰落和隐藏终端场景下单节点的检测概率惨到没法看。后来才发现问题不是能量检测本身不行而是单节点能看到的信道信息太有限。于是协作频谱感知成了绕不开的方向而集中式数据融合又是协作感知里最直接、最容易上手的架构。1.1 单节点感知到底输在哪里单个次级用户在判断主用户是否占用频段时手里只有自己接收天线上的那一段信号。如果主用户和这个用户之间恰好在深衰落或者中间隔了一栋楼接收信噪比可能低到-15dB以下这时候任何基于能量的检测器都会变成“盲猜”。更麻烦的是噪声功率本身不是一成不变的温度变化、射频前端增益漂移都会让噪声底噪上下浮动。单节点没有多余信息来消除这种不确定性所以误检和漏检很难同时压住。协作的本质是把分布在空间不同位置的多个次级用户的观测结果放到一起利用空间分集来对抗单点的深度衰落和阴影。换句话说A节点听不见的频段B节点可能听得清清楚楚B节点被干扰带偏的时候C节点还能拉一把。只要融合策略不是太离谱整体检测概率一定比最差节点好甚至能逼近最优节点的水平。1.2 集中式数据融合的两种玩法集中式数据融合里通常分硬融合和软融合两类。硬融合是每个节点先自己做出“有信号/无信号”的二元判决再把0/1结果发给融合中心最终用AND、OR或者多数表决来定案。好处是通信开销极小缺点是本地判决一旦出错融合中心很难纠回来因为信息被量化成了一比特。软融合则是节点把本地统计量直接发到融合中心比如接收能量、匹配滤波输出、循环谱特征等融合中心再按照等增益合并、最大比合并或者某种学习策略做最终判决。软融合保留了更多信息理论上性能更好代价是每节点需要回传的比特数更多。PR指数检测器在集中式融合里非常适合扮演“软统计量生成器”的角色因为算出来是一个标量既不像原始I/Q数据那样携带笨重也没有把信息压缩到只剩一比特那么极端。1.3 PR指数检测器在这个框架里扮演什么角色PR指数检测器做的事情是把每个节点接收信号的能量序列拿去做一次“分布形状体检”输出一个代表当前能量分布与噪声参考分布差异的标量。这个标量丢给融合中心融合中心把所有节点标量做加权平均或加权求和再与门限比较。整个过程很干净也不用在节点间同步复杂信息。集中式协作感知的核心诉求就是“各节点上报有效信息、中心统一决策”PR指数检测器正好卡在这个需求点上。我在实现过程中比较喜欢这种设计的另一个原因是它天然避开了对主用户信号具体调制方式的依赖。你不需要知道主用户是OFDM还是单载波甚至不需要知道信号带宽的精确边界只要能量分布的形状因为信号出现而发生变化PR指数就能反映出来。这点在工程落地时非常实用。2. Pietra-Ricci指数检测器的数学直觉与实现原理很多第一次看到“Pietra-Ricci指数”这个名字的人会以为是什么深奥的现代统计工具其实它的源头可以追溯到经济学里的不平等度量。后来做模式识别和信号检测的人发现这类指数完全可以用来衡量两个概率分布之间的分离程度于是就被搬进了频谱感知的领域。2.1 从不平等指数到分布距离Pietra指数最早用来描述收入分配的不平等程度数值越大表示资源分布越不平均。Ricci指数也类似它们实际上都跟洛伦兹曲线、平均绝对偏差这些概念有联系。但在检测问题里我们更关心的是它作为一种分布距离的性质。不同文献给出的定义虽然写法不同但在离散化实现时我最常用也最推荐的是基于概率密度L1距离的等价形式PR (1/2) * ∫ | p_obs(x) - p_ref(x) | dx其中p_obs是观测能量序列的概率密度函数p_ref是纯噪声条件下的参考概率密度函数。这个式子算出来在0到1之间两个分布完全一样就是0完全不重叠就是1。直观理解就是把两条概率密度曲线画出来算出它们之间“不重合面积的一半”面积越大说明观测到的能量分布和噪声分布差得越远也就越有把握说主用户信号来了。在实际Matlab实现里我们不会真的去做连续积分而是把能量序列丢进直方图得到离散的概率密度估计然后求所有分箱上绝对差的和再乘以箱宽。这和式子的离散版本完全对应代码写起来也不绕。2.2 把PR指数变成判决统计量现在的问题是如何用PR指数构造频谱感知的判决统计量。最直接的做法是每个节点把接收信号分成若干帧每一帧算一次能量得到一条长度为N的能量序列。比如采样率1MHz每帧取64个采样点总帧数100那么一次感知需要6400个采样点对应6.4毫秒的感知时间。这个时间在大多数认知无线电场景里都是可以接受的。有了这100个能量值就能估计观测分布。同时在系统初始化阶段先采集一段“确定没有主用户”的纯噪声数据按同样的方式计算噪声能量直方图作为参考分布。接下来对当前观测能量直方图与参考直方图计算PR指数得到本地统计量。统计量越大越倾向于“主用户存在”的假设。每个节点都把这个标量上报给融合中心融合中心再对多个节点的PR值做等增益平均、按信噪比加权或投票最后和门限比大小。这里有个很重要的点PR指数不是直接对原始信号幅值计算而是对“能量序列”的分布计算。因为要对抗噪声和信道衰落能量域的统计特性更稳定。而且分帧能量天然是正值直方图的分辨率也容易控制不会出现原始信号双峰难建模的问题。2.3 PR检测器和传统检测器差在哪很多刚接触的人会问跟能量检测、循环平稳检测比PR指数到底强在哪我用一个表格把这几个常见检测器的特点列出来方便直观对比。检测器需要噪声功率需要信号特征计算复杂度低信噪比表现主要局限能量检测需要精确估计不需要很低一般噪声不确定性敏感匹配滤波需要需要导频/前导码很低很好需要主用户先验信息循环平稳检测不需要需要循环频率高较好计算量偏大采样要求高KS检验检测器需要不需要较低较好对局部差异敏感容易过拟合PR指数检测器需要不需要较低较好直方图参数需要调节PR指数和KS检验类似都属于非参数思路但KS检验用的是经验CDF的最大差值只盯着一个点一旦直方图采样不好判决稳定性会受影响。PR指数把分布之间的差异做了积分平均平滑性更好对直方图分箱的敏感度相对低一点。当然这只是相对而言分箱仍然是个值得认真调参的地方后面我会专门讲。3. Matlab代码实现从零搭一套集中式PR检测仿真代码部分我尽量给得完整但也不是把所有工具函数都堆出来而是把核心逻辑拆开讲清楚。我自己跑仿真的时候喜欢把参数、子函数、主循环分开写这样后期改节点数、改信噪比都方便。3.1 仿真参数怎么定先定一组不容易出错的参数fs 1e6; % 采样率 1MHz M 64; % 每帧采样点数 N 100; % 帧数也是一次感知的能量样本数 P 8; % 协作节点数 numTrials 1000; % 蒙特卡洛次数 PfaTarget 0.1; % 目标虚警概率 snrDb -10; % 单节点信噪比为什么要用每帧64点因为64点做能量求和后能量序列的分布已经比较接近高斯或卡方分布直方图不会出现太多空箱。帧数100是为了平衡感知时间和分布估计精度。如果N太小比如只有20帧直方图会特别毛糙PR指数波动很大门限很难标定。如果N太大比如1000帧性能确实会更稳但感知时延可能超过认知协议的规定。8个节点是我用来演示集中式增益的典型值节点太少看不出协作优势太多则会增加回传开销和融合中心的计算压力。3.2 核心子函数计算PR指数下面这个函数是整个仿真的心脏。它接收当前节点的能量序列、参考噪声直方图和直方图边界返回一个PR值。function PR pr_detector(y, noiseHist, binEdges) y y(:); % 防止观测能量超出参考直方图边界 y(y binEdges(1)) binEdges(1); y(y binEdges(end)) binEdges(end); obsHist histcounts(y, binEdges, Normalization, pdf); delta binEdges(2) - binEdges(1); PR 0.5 * sum(abs(obsHist - noiseHist)) * delta; end注意我把超出边界的能量值强制截断到边界上这是为了避免histcounts把越界样本丢到统计范围之外。实际操作中如果参考噪声直方图的最大边界只有噪声能量最大值的1.5倍一旦主用户信号比较强观测能量很容易超过边界截断是必须的。参考直方图的生成方式如下先产生一段纯噪声能量序列再统计概率密度noiseRef zeros(N, 1); for k 1:N noise (randn(M,1) 1j*randn(M,1)) / sqrt(2); noiseRef(k) sum(abs(noise).^2); end binEdges linspace(0, max(noiseRef)*2, 51); noiseHist histcounts(noiseRef, binEdges, Normalization, pdf); noiseHist noiseHist(:);这里用复高斯噪声方差归一化到1每帧64点意味着一个能量样本的均值是64。分箱数选50箱宽大概是(max(noiseRef)*2)/50。箱子太少分布细节被抹掉箱子太多每个箱子里平均样本数太少统计波动大。50箱对100个样本来说是一个比较稳的起点。3.3 蒙特卡洛主循环与阈值标定集中式融合需要一个门限标准做法是先跑“无主用户”假设下的融合结果取虚警概率对应的分位数作为门限。我习惯先单独生成numTrials次纯噪声场景得到融合统计量的分布再用quantile求门限。TfH0 zeros(numTrials, 1); for tr 1:numTrials localPR zeros(P, 1); for p 1:P y zeros(N, 1); for k 1:N noise (randn(M,1) 1j*randn(M,1)) / sqrt(2); y(k) sum(abs(noise).^2); end localPR(p) pr_detector(y, noiseHist, binEdges); end TfH0(tr) mean(localPR); end threshold quantile(TfH0, 1 - PfaTarget);融合规则先用最简单的等增益平均也就是把每个节点的PR值求平均。等增益的好处是不需要估计信道信噪比实现在线计算量很小。当然如果各节点信噪比差异很大最好换成按SNR加权后面会提到。有信号场景的生成需要加入一条主用户信号我这里用随机复高斯信号近似一个宽带主用户按目标SNR设置幅度signal (randn(M,1) 1j*randn(M,1)) / sqrt(2); signal signal / sqrt(mean(abs(signal).^2)); noisePower 1; amp sqrt(noisePower * 10^(snrDb/10)); rx amp * signal (randn(M,1) 1j*randn(M,1)) / sqrt(2);然后同样分帧算能量生成当前节点的能量序列再用pr_detector得本地PR。最后融合中心做等增益平均和门限比较统计检测概率。3.4 融合中心决策与性能统计主循环里记录命中次数最后除以总试验次数就是检测概率。我通常把H0和H1的蒙特卡洛循环分开写避免重复生成数据导致混乱。跑完之后用roc类型的分析思路扫描一系列门限得到完整的ROC曲线而不是只给一个工作点。因为论文或者项目汇报里只有单一Pd/Pfa点说服力不够领导或者导师一定会问“曲线呢”。在Matlab里可以把threshold设为TfH0的分位数序列比如从0.01到0.99算出一组Pfa和一组Pd然后画曲线。这个操作不复杂就是多跑几轮循环但信息量比单点大得多。4. 仿真结果怎么看性能提升与坑位提醒代码跑通只是第一步真正折磨人的是结果分析和参数调优。我第一次跑出来的结果单节点检测概率只有0.3融合之后到了0.8心里挺爽。但后来改了几个参数后融合性能反而下降才意识到这里面坑不少。4.1 单节点与集中式融合的ROC对比按我上面那组参数在单节点信噪比-10dB、目标虚警0.1的条件下单节点PR检测器的检测概率大概在0.35到0.45之间具体取决于每次蒙特卡洛的噪声种子。8个节点等增益融合后检测概率能提高到0.75到0.85。这个提升幅度符合预期因为8个独立节点的统计量平均之后方差被压低分布重叠区域明显变小ROC曲线会整体向左上方移动。如果信噪比继续降到-15dB单节点基本躺平检测概率可能只有0.18。8节点融合还能勉强跑到0.5左右再想往上就得增加节点数、增加帧数或者改用按信噪比加权的融合规则。所以PR指数检测器不是银弹它靠的是协作带来的空间分集收益。4.2 不同融合规则、节点数对检测概率的影响我对比过三种融合规则等增益平均、按信噪比加权、以及硬判决OR融合。在节点间信噪比相同的情况下等增益和加权效果几乎一样OR融合会差一些因为本地判决一旦硬性阈值化部分节点在低信噪比下输出0/1时已经丢失了大量软信息。节点数从1涨到4、8、16时检测概率的增速是先快后慢。4节点相比单节点提升最明显8节点还有明显收益16节点收益就开始平缓。而且节点数越多门限标定的蒙特卡洛次数要求也越高否则门限本身的波动会影响性能比较。我建议至少跑numTrials 2000以上我有时候甚至跑5000次虽然耗时一些但曲线才光滑。4.3 阈值标定和大数定律陷阱门限标定最容易被忽略的是H0门限分布和H1检测分布必须在完全相同的参数下生成包括参考直方图、分帧方式、节点数、融合规则。如果你H0用了100帧H1用了200帧那门限直接失效。还有一点融合中心收到的本地PR值实际不是严格独立同分布的因为所有节点共用同一参考噪声直方图如果参考直方图本身有偏所有节点的统计量会带上系统性偏差。在蒙特卡洛里这不算大问题但实机测试时不同节点的射频前端噪声功率不一样必须每个节点单独建立参考分布或者做归一化处理。另外quantile在样本量不足时给出的门限很不稳定。我踩过最明显的一次是用200次蒙特卡洛去标定门限换了随机种子以后虚警概率从0.1直接跳到0.18。后来把H0试验次数加到2000门限才稳定下来。5. 工程落地中的常见问题与排查技巧从Matlab仿真到实际系统中间隔着一整条“噪声参考漂移、实时直方图、定点化”的鸿沟。这里分享几个我自己调试时经常检查的点。5.1 噪声参考不准PR指数直接漂移PR检测器依赖参考噪声分布等于把能量检测对噪声功率的依赖转移成了对噪声分布形状的依赖。如果实际噪声功率比参考分布高3dB那么观测能量序列会整体右移PR值会显著变大虚警率跟着飙升。解决思路不是不做参考而是尽量做实时参考校准。比如在认知系统认定的空白频段里定期采集噪声样本更新参考直方图。如果实在没有纯净噪声时段可以退一步用参数化的方法假设噪声能量服从Gamma分布或卡方分布根据历史样本实时估计均值和方差再生成理论参考直方图。这样虽然牺牲了一点非参数优势但至少能抗住噪声功率的慢变。5.2 帧数和直方图分箱怎么配合分箱数量是PR检测器绕不开的超参数。我的一般经验是让每个箱子平均不少于2到3个样本。如果100个能量样本分50箱是极限再多空箱就多了PR值会虚高。另外一个技巧是适当把直方图边界放宽到参考噪声最大值的2倍这样有信号时能量序列不会大量被截断到最后一个箱子尤其是强信号场景下截断会让PR值出现饱和反而掩盖不同信噪比之间的差异。分箱宽度也可以用公式来定比如K ceil(sqrt(N))这种经验法则。但我觉得最实际的做法是固定箱数跑一遍ROC看曲线变化。同一个NK30、50、80各跑一次选ROC曲线下面最大的那个。我自己最终常用N100、K50这套组合。5.3 从Matlab走向实时系统的几个改动Matlab里跑历史样本很容易但实时系统里直方图更新是个麻烦事。节点每帧算一个能量然后用滑动窗口维护最近N帧能量再实时更新直方图。Matlab的histcounts每次调用开销不小在实时处理器上建议预先计算每个能量值落在哪个箱子然后只更新对应箱子的计数PR值用增量方式更新避免每帧全量重算。另外如果融合中心要按SNR加权节点就得顺带估计本地SNR。最省事的办法是同时上报本地PR值和本地平均能量中心根据平均能量估算SNR权重。但要注意本地平均能量在有信号时也会变大直接用这个能量做权重会把强节点抬得过高不一定是最优的。可以改用长期平均的背景噪声功率来归一化。最后再分享一个我实际调试时的小技巧先在极低信噪比下跑少量试验把每节点的PR值打印出来看看是否随着SNR有单调变化。如果你发现PR值在无信号和有信号时差别很小大概率不是检测器的问题而是参考直方图生成错了比如参考能量没有归一化或者分帧方式前后不一致。先解决这个再谈融合算法能省下不少时间。