
1. 为什么单节点频谱感知不可靠从隐藏终端说起协作频谱感知这几年在认知无线电里被提的最多的一句话就是“单节点不可靠”。信号衰落、路径损耗、障碍物遮挡任何一个因素都可能让某个感知节点对主用户的存在完全“失明”。你可能会想直接在接收端量一下接收功率不就行了吗确实可以这就是经典的能量检测但在衰落信道下单个节点的瞬时信噪比可能低到完全没法判断比如主用户和感知节点之间隔了一栋楼检测概率可能掉到0.3以下。这种场景一旦发生盲目认为“频段空闲”然后让次用户接入就会对主用户产生严重干扰。解决思路其实很朴素多找几个节点一起感知把各自的观测结果提交到融合中心由融合中心综合所有信息给出一个全局判决。这就是协作频谱感知的基本思想。但这带出一个新问题——每个节点到底往融合中心传什么传原始采样波形带宽开销大到不现实只传本地0/1判决又会丢掉很多统计信息融合增益有限。于是工程上就需要一种折中方案节点本地提取一个足够反映“信号存在与否”的统计量把统计量传给融合中心再做集中式融合。Pietra-Ricci指数检测器就是这种思路下的一个典型代表。它基于接收信号的样本协方差矩阵特征值分布设计统计量天然不依赖噪声功率先验非常适合抗噪不确定性环境中的协作检测。这篇文章我打算从原理讲到Matlab实现重点拆开数据融合的门限设计和工程上容易踩的坑适合正在做协作频谱感知方向的研究生或者想快速用Matlab搭一套检测仿真的工程师。如果你时间紧可以直接跳到第4节抄代码跑数据如果你想理解为什么PR统计量能抵抗噪声不确定性、为什么融合门限必须在融合端重新标定那还是把前几节一起看完更划算。1.1 单节点检测的困境单节点检测的问题不只是“信号弱”。更麻烦的是噪声不确定性。协方差、信道、温度、器件老化都会让接收机的噪声功率随时间变化能量检测器的门限又是直接依赖噪声功率的一旦实际噪声功率和预设值偏差1dB虚警概率就可能从0.1跳到0.3。在认知无线电场景里虚警概率升高意味着频谱资源被白白浪费漏检则意味着对主用户的干扰失控两边都不能容忍。所以我个人一直觉得协作感知的首要目标不是把检测概率堆到99%而是让整个系统在噪声不确定性和信道衰落面前依然行为稳定。单个节点再怎么优化面对深度衰落和隐藏终端都是脆弱的。隐藏终端这个场景稍微展开说一下主用户发射机附近明明有人在使用频谱但因为遮挡物或阴影衰落某个感知节点完全感知不到信号。这时候如果系统只依赖这个节点做判断就会误判频段空闲次用户的接入就会干扰主用户的接收机。多节点协作的本质就是利用空间分集对冲这种单点失明。1.2 协作感知的系统模型集中式协作频谱感知的系统模型可以抽象成这样一个主用户发射机M个感知节点一个融合中心。每个感知节点独立接收信号在本地提取统计量然后通过报告信道把统计量发送给融合中心。融合中心把所有节点的统计量按某种规则合并得到一个全局判决统计量再用门限判决输出“主用户存在”或“主用户不存在”。这个模型有个很关键的特征感知信道和报告信道是分开的。很多初学者仿真的时只关心感知信道的信噪比完全忽略了报告信道噪声实际工程里报告链路本身的可靠性也是要算进系统预算里的。第6节我会专门讲这部分坑。在这之前我们先默认报告信道是理想的先把核心检测和融合逻辑跑通。2. Pietra-Ricci指数检测器特征值离散度如何暴露信号如果只用一句话概括PR指数检测器的核心那就是噪声环境下的协方差矩阵特征值趋于均匀混入相关信号之后特征值会变得不均匀PR指数就是用来量化这种不均匀程度的。这个思路最早并不在通信领域而是统计学里用来衡量收入分配不平等程度的指标后来被信号处理研究者借过来用在频谱感知里意外地好用。2.1 从接收信号到样本协方差矩阵考虑一个感知节点在某段时间内采集的复基带信号。如果节点有多根天线每一次采样天然形成一个多通道向量如果节点只有单天线也可以通过滑窗或者过采样构造“虚拟多通道”向量。这里假设每个观测向量是L维的一共N个快拍组成一个L行N列的矩阵X。样本协方差矩阵就是R X * X / N其中X表示共轭转置。这个R是L乘L的方阵数学上表征了各通道之间的相关结构。在纯噪声环境下各通道之间独立R在主对角线附近比较均匀特征值大致相等。如果有主用户信号混进来并且信号在通道之间存在时间或空间相关性R的特征值就会拉开差距——少数大特征值对应信号子空间其余小特征值对应噪声子空间。我实际做仿真时最常用的构造方式是对单天线接收序列做不重叠滑窗。把长度为L*N的接收序列y排成L行N列的矩阵让矩阵第n列恰好是原始序列的第((n-1)L1)到(nL)个采样点。这样相邻采样点被安排到同一列的不同“虚拟通道”里AR信号或窄带信号在相邻采样点间有相关性而白噪声没有PR检测器就能把这两者区分开。这个构造方式非常关键第6节我会专门说为什么要这么做以及怎么做容易出错。2.2 PR指数的定义与性质有了协方差矩阵R接下来计算它的L个特征值按从大到小排列。PR指数的工程实现形式我常用下面这个T_PR mean(lambda.^2) / mean(lambda)^2 - 1展开写就是T_PR等于“特征值平方的均值”除以“特征值均值的平方”再减1。这个量有两个很直观的性质。第一如果所有特征值完全相等T_PR等于0特征值差异越大T_PR越大。第二它只取决于特征值之间的相对比例如果特征值整体乘以一个常数c比如噪声功率变大分子分母同时乘c²比值不变。这就是PR检测器抗噪声不确定性的数学本底。判决规则很简单T_PR超过门限gamma判定主用户存在也就是H1否则判定H0。门限怎么定第3节和第4节都会涉及。需要强调的是T_PR这个值本身不是概率也不是似然比它只是一个特征值离散度的度量所以门限必须在给定虚警概率下重新标定绝对不能拿一个经验常数到处套。2.3 PR检测器、能量检测器与MME检测器的定位差异做检测器性能对比时最常被拉出来比较的就是能量检测ED、最大最小特征值比检测MME和PR指数检测。能量检测的道理最简单看接收信号的总功率。它的判决统计量依赖噪声功率门限需要知道当前噪声到底有多大这在噪声不确定时会出问题。MME则是看最大特征值和最小特征值的比值不需要知道噪声功率但对小特征值的估计误差很敏感当L较大而N不够大时最小特征值估计不准性能波动明显。PR指数介于两者之间既不要噪声功率先验又不像MME那样只盯着最大和最小两个极端值而是把全部L个特征值都纳入统计。它对于“多个特征值被抬高但最大特征值不明显”的情况更鲁棒。我用一句话总结这三个检测器的差异给刚接触这个方向的朋友参考ED在理想环境下最好MME在抗最坏情况上很激进但N不够时容易翻车PR在稳健性和计算复杂度之间比较平衡。3. 集中式融合多节点PR统计量怎么合成全局判决本地每个节点都算出一个PR统计量之后问题就变成融合中心拿到M个T值怎么合成一个全局判决这是整个系统设计里最容易被低估的一步因为融合规则直接决定性能上限而这个上限和检测器的选择是正交的。检测器再强融合规则不对增益照样发挥不出来。3.1 软融合与硬融合的取舍融合规则大体分两类软融合和硬融合。硬融合是每个节点先把本地统计量和一个本地门限比较产生0/1判决再把这1bit上传给融合中心融合中心按“至少K个节点判有”或者“所有节点都判有”之类的规则输出。优点是报告信道开销极小、抗报告链路噪声能力强缺点是本地判决已经把统计量里的幅度信息扔掉了融合中心没法拿不同节点的统计量置信度做加权。软融合则更“奢侈”一些每个节点上传本地的连续统计量T_m融合中心做加权或直接平均再做一次全局门限判决。这样做保留了更多信息性能上限更高但对报告信道质量敏感。标题既然写的是集中式数据融合我默认走的就是软融合路线这也是研究论文里更常见、能拉开性能差距的场景。我给的Matlab代码采用等增益合并也就是把M个节点的T_m直接取平均因为实现最简单、参数最少用来验证基本逻辑最合适。如果你想做加权融合一种很自然的扩展是按噪声功率反比加权信噪比高的节点权重高信噪比低的节点权重低。但这样做要求每个节点能估计本地噪声功率在噪声不确定性场景下引入的误差可能比权重带来的增益还大。工程上我见过有人直接忽略噪声功率估计单纯用等增益合并性能也挺稳。所以我的建议是第一版系统别急着上复杂的权重方案等增益合并跑通后再考虑要不要加报告信道噪声模型和加权策略。3.2 门限定标必须在融合之后初学者经常犯的一个错误是先在本地把每个节点的门限算好然后把本地判决结果融合。或者更隐蔽一点先对单个节点的T_PR做门限判决假设融合中心只是“多数表决”这是硬融合的思路和软融合系统根本不是一回事。在软融合的框架下全局门限必须基于“融合之后”的统计量在H0条件下的分布来确定。为什么因为即使每个节点的T_PR在H0下的分布都相同M个独立统计量的均值分布和单个统计量的分布也不是同一个东西。均值更加紧缩方差小得多直接沿用单节点门限会导致实际虚警概率远低于目标值看起来好像“性能很好”其实是把门限定松了检测概率也被压低了虚警和检测都不在标称点上。我自己调试的时候遇到过类似情况当时第一反应是“融合增益这么大吗”后来检查门限才发现是标定位置错了虚警概率只有0.01而不是目标0.1属于典型的虚假性能提升。正确做法是在H0条件下跑出融合后统计量T_fus的大量样本取置信水平1-Pf对应的分位数作为门限。这个流程看起来繁琐但做一次就行第4节代码里就是按照这个思路实现的。4. Matlab实现从零搭起PR协作检测仿真这一节把完整实现流程拆开讲。代码我尽量写得直白方便直接改成你自己项目里的模块。需要说明的是下面的仿真中我用了不重叠滑窗构造虚拟通道这是单天线节点构造样本协方差矩阵最常见的方式。4.1 仿真参数与系统模型先定义系统参数。我这里设M4个协作节点每个节点的滑窗长度L4个采样点快拍数N128也就是说每个节点实际采集L*N512个采样点。主用户信号用复高斯白噪声经过一阶AR滤波器产生目的是让信号在时间上具备一定的相关性——这是PR检测器能工作的前提之一。每个节点的信道是独立瑞利衰落噪声是复高斯白噪声感知信噪比从-15dB到0dB扫描。clear; clc; close all; M 4; % 协作节点数 L 4; % 滑窗长度等价于虚拟通道数 N 128; % 快拍数构造L*N的观测矩阵 snr_vec -15:3:0; % 各节点平均感知SNR(dB) Pf_target 0.1; % 目标虚警概率 nMc 2000; % 蒙特卡洛次数 rng(2025); % 固定随机数种子便于复现信号模型这块主用户信号s的长度是L*N经过AR(1)滤波后再归一化到单位功率。这样设置的好处是不同节点的噪声方差可以按SNR公式直接算不用额外校准信号功率。4.2 核心函数PR统计量的实现PR统计量封装成一个独立的Matlab函数输入任意L行N列的复矩阵输出一个标量T。这个函数内部就是协方差矩阵、特征值分解、二阶矩归一化三步没有任何多余逻辑。建议你在自己的项目里把这个函数单独存成一个文件因为后面无论是单个节点检测还是融合后检测都要反复调用它。function T prDetector(X) % X: L x N 复基带观测矩阵 % 输出: PR统计量 T mean(lambda.^2)/mean(lambda)^2 - 1 L size(X, 1); N size(X, 2); R (X * X) / N; % 样本协方差矩阵LxL lambda eig(R); % 特征值按升序返回但PR只依赖整体分布 lambda lambda(:); mu1 mean(lambda); % 特征值均值 mu2 mean(lambda.^2); % 特征值平方均值 T mu2 / (mu1^2) - 1; end这里有个小细节值得注意eig函数默认返回的特征值顺序不保证从大到小但PR统计量只用到全体特征值的均值和平方均值顺序不影响结果所以不需要sort。如果哪位同学是在这个函数基础上做MME检测对比那就要改成求最大特征值和最小特征值到时候再排序也不迟。4.3 门限生成与蒙特卡洛检测循环门限生成思路就是纯噪声下反复调用prDetector拿到融合统计量T_fus的大量样本然后取分位数。纯噪声样本不乘任何噪声系数因为PR统计量天然与绝对噪声功率无关仿真里直接生成单位方差复高斯噪声就行少一个变量简化不少。% 在H0(纯噪声)下生成融合统计量的分布 T0 zeros(1, nMc); for k 1:nMc noise (randn(L, N) 1i*randn(L, N)) / sqrt(2); T0(k) prDetector(noise); end gamma quantile(T0, 1 - Pf_target); % 目标虚警概率对应的门限检测概率循环部分每个蒙特卡洛实验里先生成一段主用户信号然后每个节点独立乘衰落系数、叠加对应SNR的噪声计算本地PR统计量最后在融合中心取平均得到T_fus再与门限比较。Pd zeros(size(snr_vec)); for idx 1:length(snr_vec) snr snr_vec(idx); sigma_n sqrt(10^(-snr/10)); % 由SNR反推噪声标准差 cnt 0; for mc 1:nMc % 生成主用户信号AR(1)过程单位功率 s 0.5 * (randn(1, L*N) 1i*randn(1, L*N)); s filter(1, [1 -0.9], s); s s / sqrt(mean(abs(s).^2)); T_fus 0; for m 1:M h (randn 1i*randn) / sqrt(2); % 节点独立瑞利衰落 noise sigma_n * (randn(1, L*N) 1i*randn(1, L*N)) / sqrt(2); y h * s noise; % 接收信号 X reshape(y, L, N); % 滑窗构造虚拟通道 T_local prDetector(X); T_fus T_fus T_local; end T_fus T_fus / M; % 等增益融合 if T_fus gamma cnt cnt 1; end end Pd(idx) cnt / nMc; end这段循环的结构分得很清楚信号生成、每节点感知、融合判决三层嵌套。实际跑的时候nMc2000、M4、SNR扫描7个点大约几分钟能跑完。如果想快速验证代码逻辑可以把nMc先降到500趋势不会变就是曲线毛糙一点。4.4 绘制ROC与Pd-SNR曲线检测器仿真除了固定Pf点看Pd随SNR变化更完整的做法是画出ROC曲线也就是固定SNR扫不同的门限得到一组Pd-Pf点。门限扫描不能瞎编还是同样的思路——H0下生成T_fus的所有样本然后依次取不同分位数作为门限。% ROC曲线固定SNR-10dB扫门限分位数 snr_fix -10; sigma_n sqrt(10^(-snr_fix/10)); T0_roc zeros(1, nMc); for k 1:nMc noise (randn(L, N) 1i*randn(L, N)) / sqrt(2); X reshape(noise, L, N); T0_roc(k) prDetector(X); end gammas quantile(T0_roc, linspace(0.01, 0.99, 50)); Pf_roc zeros(size(gammas)); Pd_roc zeros(size(gammas)); for gi 1:length(gammas) gamma_now gammas(gi); pf_cnt 0; pd_cnt 0; for mc 1:nMc % H0样本 noise (randn(L, N) 1i*randn(L, N)) / sqrt(2); X0 reshape(noise, L, N); T0_now mean(arrayfun((m) prDetector(X0), 1:M)); % 简化示意 if T0_now gamma_now pf_cnt pf_cnt 1; end % H1样本 s 0.5 * (randn(1, L*N) 1i*randn(1, L*N)); s filter(1, [1 -0.9], s); s s / sqrt(mean(abs(s).^2)); T1_now 0; for m 1:M h (randn 1i*randn) / sqrt(2); noise_m sigma_n * (randn(1, L*N) 1i*randn(1, L*N)) / sqrt(2); X1 reshape(h * s noise_m, L, N); T1_now T1_now prDetector(X1); end T1_now T1_now / M; if T1_now gamma_now pd_cnt pd_cnt 1; end end Pf_roc(gi) pf_cnt / nMc; Pd_roc(gi) pd_cnt / nMc; end上面这段ROC代码里我用了arrayfun简化写法实际跑的时候为了效率建议用循环替代。ROC曲线跑起来明显比单点Pd-SNR要慢因为每个门限点都要重新蒙特卡洛。工程上常用的做法是先只跑几个关键门限点比如Pf0.01、0.05、0.1、0.2连成折线足够说明趋势。5. 仿真结果与关键参数的影响规律我自己按上面代码跑的典型结果先说明数值会有随机波动但趋势是稳定的。在Pf_target0.1、M4、L4、N128、nMc2000的参数组合下单节点SNR-15dB时检测概率大约只有0.12说明这个信噪比已经接近检测极限SNR-10dB时Pd可以到0.55左右SNR-7dB时超过0.85到SNR-3dB以上基本就是0.98以上了。加了四个节点融合后相比单节点同一个SNR点检测概率普遍能提升0.2到0.4尤其在中等信噪比区域提升最明显。5.1 节点数M与滑窗长度L的取舍节点数M从1增加到8检测性能确实是单调提升的但边际增益递减。M从1到4的增益非常可观M从4到8的增益明显变缓。原因也很好理解每增加一个节点融合中心就多一个独立观测置信度自然增加但多个节点的统计量取平均之后方差随M的增大是反比关系方差从大变小带来的边际收益在M足够大后会变弱。所以工程上不是节点越多越好还要考虑报告信道带宽和系统成本通常在4到8个节点之间选。滑窗长度L这个参数比很多人想的敏感。L过小比如只有2特征值数量太少PR统计量对“特征值不均匀”的描述能力有限L过大比如取16在总采样长度不变的情况下N会等比例缩小协方差矩阵的估计质量恶化检测概率反而下降。我的经验是L取4到8之间比较稳妥。如果想调性能优先调N而不是调L因为协方差矩阵估计质量基本由N决定N越大T_PR在H0下的方差越小门限可以压得更紧H1下的统计量也能更稳定地抬起来。5.2 与能量检测器的性能对比在这个场景下我做过的对比实验也值得提一下。当噪声功率被精确知道时能量检测器的性能确实略优于PR检测器因为能量检测把接收信号的全部能量都利用了而PR检测器还要经过特征值变换丢掉了一部分能量信息。但问题是“精确知道噪声功率”在真实系统里几乎不存在。我做过一个简单的实验仿真时把噪声功率估计值故意偏移±1dB能量检测器的曲线立刻变得很难看在SNR-12dB以下几乎失去意义PR检测器基本不受影响因为在特征值比值计算中噪声功率这一项被消掉了。如果你以后要写工程报告或者论文我建议把噪声不确定性这个维度加进对比实验这比单纯比较理想条件下的Pd更能体现PR检测器的真正价值。很多论文里说能量检测器“抗噪声不确定性差”其实只有在实验里让噪声功率偏移才能把这个结论落到数据上否则读者感受不到差异。5.3 门限分位数与虚警概率的联动还有一个容易忽略的点是门限gamma取的是H0分布的分位数Pf_target设得越小门限越大检测概率也会跟着下降。这是检测器内部不可避免的trade-off。如果你发现实验结果里Pd始终上不去先去看看Pf_target设的多少有时候不是检测器不行而是你把虚警概率压得太低了。认知无线电系统通常会根据主用户保护的严格程度设定Pf一般在0.05到0.2之间不必追求越低压越好。6. 工程实现中的坑位与绕坑方法这套仿真代码看起来不长但我在实际调试时踩过的坑不比写代码的时间少。挑几个最容易影响结果的问题集中说一下算是给后来者排雷。6.1 协方差矩阵的秩与正定性问题样本协方差矩阵R是L乘L的矩阵要保证它满秩一个基本条件是N不小于L。但这只是理论上限实际使用里N必须明显大于L才够稳。当N和L接近时小特征值会被严重低估导致T_PR在H0下方差变大、门限抖动最直接的后果就是ROC曲线不稳定。如果你跑出来的曲线一次一个样先检查是不是N/L太小。我的默认建议是N至少是L的10倍以上N128对应L不超过8这也是我上面参数选L4而不是L16的原因之一。6.2 reshape滑窗构造虚拟通道的顺序别搞反Matlab的reshape函数是按列填充的也就是说对于一个长度为L*N的行向量yreshape(y, L, N)生成的矩阵第一列是y的前L个元素第二列是接下来的L个元素。如果你心里想的是“第一行是连续N个时间样本”这种构造方式会完全错位导致所谓的“虚拟通道”之间根本没有你想要的相邻采样相关结构检测性能自然崩掉。我自己第一次实现时就栽在这里当时还以为是自己公式写错了后来手推了一个小例子才发现是reshape顺序问题。要验证构造方式对不对很简单生成一个线性序列1到L*Nreshape成矩阵后打印第一列看它是不是[1;2;...;L]。如果是说明相邻采样被安排到了同一列的不同行这才是我们想要的虚拟通道结构。6.3 报告链路噪声与量化误差的干扰前面所有仿真都假设报告信道理想节点往融合中心传的统计量分毫不差。真实系统里这是不可能的。PR统计量本身是个连续值要想数字传输必然要做量化。量化bit数太少比如只保留1bit那其实等价于硬融合软融合的性能优势就全丢了量化bit数太多报告链路开销又变大。我实践下来的经验是软融合最少也要3bit量化最好5bit以上否则统计量的幅度信息损失严重。报告链路本身也可能加噪声。如果报告链路的SNR比较高比如高于10dB对检测性能影响很小如果报告链路质量差还不如改用硬融合更划算。这一点在系统设计时经常被忽略建议你在仿真后期加上一个简单的报告链路模型在T_m上叠加一个零均值高斯扰动扰动方差按报告链路的SNR反算然后重新做融合、重新标定门限。整个过程跑完之后你会发现门限又要重算一遍——这就对了任何链路模型改变后门限都必须跟着重标定这是协作频谱感知仿真里最需要刻进肌肉记忆的一条规则。我做这套仿真最深的体会是PR检测器本身的实现并不复杂真正费时间的往往是你没有意识到哪些参数和环境变了之后门限必须重算。只要把这个“改参数就重标门限”的习惯养好后面不管是换信号模型、换信道模型还是加报告链路都只是往这套框架里填模块的事情。希望这篇内容能帮你少熬几个晚上如果你跑出来的趋势和我不一样先别急着怀疑代码优先检查协方差矩阵构造顺序和门限标定位置八成问题就出在那里。