做无线通信仿真的人八成都被频谱感知里的能量检测坑过。它简单、直观但一遇到噪声不确定度就露馅。我最近在复现集中式协作频谱感知时发现主用户信号在多个认知节点之间是有空间相关性的融合中心只要把各节点采样拼成矩阵这种结构就藏在协方差矩阵里。Pietra-Ricci指数检测器PR检测器正好是专门干这个的它不直接估计噪声功率而是靠协方差矩阵特征值的“扁平程度”来做判决。这篇文章先讲清楚PR指数的原理再给出完整的Matlab仿真代码最后聊几句调参和落地时容易踩的坑。适合刚开始接触协作频谱感知或者想在Matlab里把特征值类检测器跑通的同学参考。1. 从能量检测到结构检测为什么集中式协作感知需要PR指数1.1 能量检测的“噪声不确定度墙”先说个最基础的问题为什么非要用PR指数很多资料一上来就讲故事但我觉得得先从能量检测的痛点讲起否则你根本不会珍惜PR检测器。能量检测的检验统计量是接收信号在感知时隙内的平均功率[ T_E \frac{1}{M}\sum_{m1}^{M}|y(m)|^2 ]判决时拿它和门限比门限通常写成 (\gamma \sigma_w^2 \cdot \lambda)其中 (\sigma_w^2) 是噪声方差。麻烦就在这噪声方差是估计出来的而且会随着温度、带宽、射频前端增益变化。实际测量中噪声功率波动 12 dB 非常正常到了低信噪比区域检测性能会急剧恶化这就是著名的“信噪比墙”问题。协作频谱感知能缓解一部分问题。多个认知用户把信息送到融合中心融合中心做判决空间分集会带来增益。可如果融合中心处理的仍然只是各节点上报的能量值那么噪声不确定度的问题只是被摊薄并没有从根本上消除。节点1的噪声估高了节点2的噪声估低了融合中心对这些估计误差其实很头疼。因此我们需要一种在理想情况下与噪声功率绝对无关的检测量。1.2 协作数据融合带来的空间相关性集中式协作频谱感知的架构非常直接每个认知用户在一个感知时隙内采样然后把采样数据或经量化的采样数据发给融合中心。融合中心拿到的是 (N) 个用户的观测序列可以拼成一个 (N \times M) 的矩阵[ \mathbf{Y} \begin{bmatrix} y_1(1) y_1(2) \cdots y_1(M) \ y_2(1) y_2(2) \cdots y_2(M) \ \vdots \vdots \ddots \vdots \ y_N(1) y_N(2) \cdots y_N(M) \end{bmatrix} ]当主用户不存在时(\mathbf{Y}) 每一行都是独立的高斯白噪声行与行之间没有任何相关性。当主用户信号出现时情况就不同了主用户信号经过不同信道到达各认知用户这些信号都来自同一个源所以不同行之间会存在相关性。你不需要精确知道信道系数是多少只要知道“相关性的存在”这件事本身就能判断主用户是否活跃。融合中心怎么量化这种相关性一个自然的选择是计算样本协方差矩阵[ \mathbf{R} \frac{1}{M}\mathbf{Y}\mathbf{Y}^H ](\mathbf{R}) 是一个 (N\times N) 的矩阵。如果只有噪声(\mathbf{R}) 应该近似为一个对角阵而且对角线上的值都接近噪声功率如果存在主用户信号(\mathbf{R}) 的非对角元素会显著增加特征值结构也会发生变化。1.3 PR指数如何抓住“结构变化”PR指数最早并不是给频谱感知用的它更多是用于衡量两个概率分布之间的差异。到了频谱感知场景里我们关心的分布是协方差矩阵特征值的分布。给定样本协方差矩阵 (\mathbf{R})其特征值为 (\lambda_1, \lambda_2, \ldots, \lambda_N)。纯噪声情况下所有特征值都挤在噪声功率附近非常“扁平”存在主用户信号时最大的几个特征值会被信号拉伸特征值分布变得“尖锐”。PR指数要做的就是把这种扁平/尖锐程度浓缩成一个标量。我在仿真中采用的PR检测统计量定义如下[ T_{PR} \frac{\mathrm{tr}^2(\mathbf{R})}{\mathrm{tr}(\mathbf{R}^2)} \frac{\left(\sum_{i1}^{N}\lambda_i\right)^2}{\sum_{i1}^{N}\lambda_i^2} ](H_0)主用户不存在下(\lambda_i \approx \sigma_w^2)所以分子约为 (N^2\sigma_w^4)分母约为 (N\sigma_w^4)因此 (T_{PR} \approx N)。(H_1)主用户存在下少数大特征值占据主导分母增长得比分子快因此 (T_{PR}) 小于 (N)。你看这个统计量里根本没有噪声功率 (\sigma_w^2)它取决于特征值之间的相对关系。这就是PR指数检测器能对抗噪声不确定度的根源。2. PR指数检测器的原理拆解与门限设计2.1 从协方差矩阵特征值看信号有无刚才用了“扁平”和“尖锐”这种比较形象的说法现在落到矩阵特征值上再讲透一点。把 (\mathbf{R}) 做特征值分解[ \mathbf{R} \mathbf{U}\boldsymbol{\Lambda}\mathbf{U}^H ](\boldsymbol{\Lambda}) 是对角矩阵对角线是特征值。在理想情况下(H_0)(\mathbf{R} \sigma_w^2 \mathbf{I}_N)特征值全部相等(\lambda_1 \lambda_2 \cdots \lambda_N \sigma_w^2)。(H_1)(\mathbf{R} \mathbf{R}_s \sigma_w^2 \mathbf{I}_N)(\mathbf{R}_s) 是信号成分的协方差矩阵。由于主用户信号占用的空间维度通常远小于 (N)所以(\mathbf{R}_s) 的秩较低叠加之后只有少数几个特征值明显变大。举个例子假设 (N4)纯噪声时特征值可能是 ([0.98, 1.02, 1.01, 0.99])都围绕着噪声方差 (1) 波动信号存在时特征值可能变成 ([4.3, 0.97, 1.03, 0.96])。这个时候计算 (T_{PR})纯噪声情况约等于 (4)信号存在时明显小于 (4)。所以PR检测器的判决逻辑可以写成[ T_{PR} \mathop{\gtrless}_{H_1}^{H_0} \gamma ]更具体地说如果 (T_{PR} \gamma)判为 (H_1)否则判为 (H_0)。2.2 检测统计量与判决准则的形式化假设我们面对的是二元假设检验问题[ \begin{cases} H_0: \mathbf{Y} \mathbf{N} \ H_1: \mathbf{Y} \mathbf{H}\mathbf{s} \mathbf{N} \end{cases} ]其中 (\mathbf{Y}) 是融合中心收到的 (N\times M) 观测矩阵(\mathbf{H}) 是信道系数向量/矩阵(\mathbf{s}) 是主用户信号向量(\mathbf{N}) 是噪声矩阵。融合中心先计算样本协方差矩阵 (\mathbf{R})再计算 (T_{PR})。由于 (T_{PR}) 是特征值的对称函数它天然不依赖特征向量所以对信道相位不敏感。这个性质很好毕竟频谱感知阶段我们通常不知道主用户信号的相位和精确信道状态信息。与经典的最大最小特征值检测器MME相比[ T_{MME} \frac{\lambda_{\max}}{\lambda_{\min}} ]MME只用了最大和最小两个特征值中间的信息全丢了。PR检测统计量用了全部特征值相当于把所有特征值的均匀度都考虑进去。在样本数不够大、特征值扰动比较明显的场景下PR检测器的统计稳定性通常优于MME。当然这个优势不是绝对的但至少从工程实现角度看PR指数的计算不需要显式求解特征值直接用矩阵的迹就能算出来数值上更稳定复杂度也更低。2.3 门限怎么定用蒙特卡洛给检测器“开门”PR检测器的理论门限推导需要用到随机矩阵理论中的Marchenko-Pastur分布或者Tracy-Widom分布推导过程相当绕而且依赖高斯噪声假设。如果你只是做工程验证或者课程仿真我建议直接用蒙特卡洛标定门限。思路很简单在 (H_0) 下独立重复生成大量噪声数据计算每次的 (T_{PR})得到 (T_{PR}) 在纯噪声下的经验分布。给定目标虚警概率 (P_{fa})取该分布的第 (P_{fa}) 分位数作为门限。为什么取左分位数因为 (H_0) 下 (T_{PR}) 分布在较大值附近(H_1) 下 (T_{PR}) 会往小值方向移动。我们要找一个门限使得纯噪声下 (T_{PR}) 小于这个门限的概率正好是 (P_{fa})。所以用quantile(T_h0, Pfa)就对了。核心代码如下numMC 5000; M 1000; N 4; Pfa 0.05; T_h0 zeros(1, numMC); for mc 1:numMC Y (randn(N, M) 1j*randn(N, M)) / sqrt(2); R (Y * Y) / M; T_h0(mc) trace(R)^2 / trace(R*R); end gamma quantile(T_h0, Pfa);这段代码看起来简单但有两个细节很容易翻车。第一(H_0) 和 (H_1) 的数据必须用独立的随机流生成不能顺手复用同一批随机数第二蒙特卡洛次数太少时门限波动很大后续仿真出来的检测概率曲线就不平滑。我一般至少用 5000 次如果目标 (P_{fa}) 是 0.01 以下建议加到 20000 次以上。2.4 PR指数与噪声不确定度的关系有人可能会问你计算 (\mathbf{R} \frac{1}{M}\mathbf{Y}\mathbf{Y}^H)里面不是也包含噪声吗怎么就说与噪声功率无关了关键在于 (T_{PR}) 的表达式。在 (H_0) 下 (\mathbf{R} \approx \sigma_w^2 \mathbf{I}_N)那么[ T_{PR} \frac{(N\sigma_w^2)^2}{N\sigma_w^4} N ]如果噪声功率从 (1) 变成 (2)分子分母同步缩放最终结果还是 (N)。也就是说门限不需要随着噪声功率调整。这个性质在实际中非常重要因为只要接收机工作在近似白噪声的频带内就不需要频繁重新标定门限。当然如果噪声是色噪声或者各节点噪声功率不一致(\mathbf{R}) 的H0结构就不再是纯比例单位阵门限会偏移。这时需要先做噪声白化预处理或者重新在相应噪声条件下标定门限。这一点留到后面“工程落地”部分再细说。3. Matlab完整实现从单次检测到蒙特卡洛曲线3.1 仿真系统参数与模型设置在写完整代码之前先把仿真参数列清楚。下面这套参数是我实际调试时用的兼顾了性能和运行时间。参数符号取值说明认知用户数N4集中式融合节点数每节点采样数M1000感知时隙内采样点数蒙特卡洛次数numMC5000用于门限标定和性能统计目标虚警概率Pfa0.05门限对应的虚警水平信噪比范围SNR_dB-20:2:0 dB仿真检测概率曲线的横轴主用户信号sQPSK单位平均功率每符号实虚部均为±1/√2噪声n复高斯白噪声单位功率实虚部独立方差各0.5信道h瑞利平坦衰落每节点一个复增益平均功率1这里有一个需要提前说清楚的信噪比定义接收端信号平均功率除以噪声平均功率。由于噪声功率归一化为 1所以给定 SNR_dB 后信号幅度就是 (\sqrt{10^{SNR_dB/10}})。如果信道是单位平均功率的瑞利衰落那么接收信号平均功率就等于 SNR。3.2 核心函数样本协方差与PR统计量为了代码整洁我把PR统计量封装成函数function T prStat(Y) % PR_Pietra-Ricci 指数检测统计量 % Y: N x M 复基带观测矩阵 % T tr(R)^2 / tr(R^2), R Y*Y/M R (Y * Y) / size(Y, 2); T trace(R)^2 / trace(R * R); end为什么不用cov(Y)因为cov默认会减去每行均值而频谱感知模型里噪声和信号通常假设零均值减均值反而会引入额外估计误差。直接用Y*Y/M更贴近理论模型也避免了样本均值不为零时的偏差。3.3 门限标定与检测概率仿真脚本下面是完整的仿真脚本框架你可以直接复制到Matlab里运行% PR_CSS_sim.m -- Pietra-Ricci指数检测器集中式协作频谱感知 clear; clc; close all; % 参数设置 N 4; % 认知用户数 M 1000; % 每用户采样点数 numMC 5000; % 蒙特卡洛次数 Pfa 0.05; % 目标虚警概率 SNR_dB -20:2:0; % 信噪比范围 % -------- 门限标定纯噪声 -------- T_h0 zeros(1, numMC); for mc 1:numMC Y (randn(N, M) 1j*randn(N, M)) / sqrt(2); T_h0(mc) prStat(Y); end gamma quantile(T_h0, Pfa); % -------- 检测概率仿真 -------- Pd zeros(size(SNR_dB)); for k 1:length(SNR_dB) snr 10^(SNR_dB(k)/10); cnt 0; for mc 1:numMC % 瑞利平坦衰落信道单位平均功率 h (randn(N, 1) 1j*randn(N, 1)) / sqrt(2); % QPSK主用户信号单位平均功率 s (sign(randn(1, M)) 1j*sign(randn(1, M))) / sqrt(2); % 噪声 noise (randn(N, M) 1j*randn(N, M)) / sqrt(2); % 接收信号 Y h * (sqrt(snr) * s) noise; % 检测 T prStat(Y); if T gamma cnt cnt 1; end end Pd(k) cnt / numMC; end % -------- 绘图 -------- figure; plot(SNR_dB, Pd, -o, LineWidth, 1.5, MarkerSize, 5); grid on; xlabel(SNR (dB)); ylabel(检测概率 P_d); title(PR指数检测器在集中式协作频谱感知中的性能);这段代码并不长但已经涵盖了门限标定、信号生成、检测判决和蒙特卡洛统计。运行一次在我的机器上大约需要几分钟如果觉得慢可以把numMC降到 2000曲线会略抖但趋势还在。3.4 怎样把仿真结果画成规范曲线在实际报告中光有单条 (P_d)-SNR 曲线不够通常要把不同 (N) 或不同 (M) 的曲线画在一起才看得出协作增益。我在仿真时会额外加一个基线单节点能量检测。单节点能量检测可以这样写% 单节点能量检测基线已知噪声方差为1 Y1 (randn(M, 1) 1j*randn(M, 1)) / sqrt(2); E sum(abs(Y1).^2) / M; % 门限需要根据噪声方差计算这里噪声方差1 gamma_energy 1 sqrt(2/M) * erfinv(1 - 2*Pfa); % 高斯近似然后把不同 (N) 值的 PR 检测器曲线叠加到同一张图上横坐标 SNR纵坐标 (P_d)。从图上你会非常直观地看到(N2) 比 (N1) 强但增益不是线性增长的(N4) 到 (N8) 的增益开始饱和。这说明协作感知并不是节点越多越好节点过多还会带来同步与回传开销需要在系统设计时权衡。3.5 集中式融合和分布式融合在代码上的区别题目里专门强调了“集中式数据融合”。刚上面的代码就是集中式融合中心直接拿到所有节点的观测矩阵 (\mathbf{Y})。如果改成分布式融合各节点只上传本地统计量而不是原始采样那么代码会有两个变化第一各节点先独立计算自己的 (T_{PR}^{(i)})对 (N1) 的协方差矩阵其实退化成能量检测第二融合中心把各节点的统计量加权合并例如[ T_{fusion} \sum_{i1}^N w_i T_{PR}^{(i)} ]由于 (N1) 时PR统计量恒等于1所以分布式场景下直接用PR指数是没意义的得用其他局部统计量如能量或循环平稳特征再在融合中心做软合并。这恰恰说明了集中式数据融合对PR检测器的重要性只有把多节点观测拼成矩阵协方差矩阵才会呈现空间结构PR统计量才能发挥作用。4. 仿真结果解读与参数调优经验4.1 不同信噪比下检测概率的典型趋势我按上面的参数跑出来的典型结果是在 (P_{fa}0.05)、(N4)、(M1000) 时SNR 从 -20 dB 增加到 -14 dB 左右时(P_d) 会从接近 0 快速拉升到接近 1。拐点越陡说明检测器性能越好。为什么会有这个拐点因为协方差矩阵估计的误差是随着 SNR 变化的。低 SNR 时信号被噪声淹没特征值结构接近纯噪声(T_{PR}) 和门限接近随着 SNR 提高信号子空间逐渐从噪声子空间“挤”出来特征值离散度增大检测概率迅速上升。如果对比单节点能量检测在同等的噪声不确定度条件下PR检测器的拐点会提前约 23 dB。这个增益来自空间协作不是算法魔法。如果你把 (N) 提高到 8拐点可能再提前 12 dB但边际增益递减。4.2 参数N和M的选择技巧参数选择是仿真最容易纠结的地方。我个人的经验是看两个比值(N) 决定空间维度(M) 决定协方差矩阵估计质量两者之间靠 (M/N) 关联。(M) 太小样本协方差矩阵偏离真实协方差矩阵(H_0) 下特征值离散度也会变大门限对应的虚警概率会高于目标值。(N) 太大需要同步的节点太多主用户信号在不同节点之间的相关性可能下降比如分布在不同位置的节点看到不同的信道衰落协方差矩阵的非对角结构不再明显。经验法则(M \ge 10N) 通常是一个比较稳的起点。如果 (M1000)(N) 最多设到 10 左右如果你需要模拟大规模节点就要相应增加采样点数。下表是我在不同参数组合下的观察NMM/N能达到的稳健程度2500250门限较稳但空间增益有限41000250推荐组合性能和开销平衡81000125需要更高信噪比才能稳定估计164000250增益饱和同步复杂度明显增加4.3 门限设置不当会怎样门限标定是整个仿真里最容易“看起来正确、实际上错误”的环节。如果你把门限设得太大比真实门限更靠右那么 (H_1) 下 (T_{PR}) 小于门限的概率变大虚警概率也随之增大。如果你把门限设得太小检测概率会下降主用户信号会被漏检这在频谱感知里更危险因为漏检意味着你会去占用一个正在被使用的频段对主用户造成干扰。蒙特卡洛标定门限时还有一点容易被忽略quantile只是经验分位数当numMC不够大时经验分位数本身的方差不可忽视。尤其目标 (P_{fa}0.01) 时如果只跑 1000 次 (H_0)那门限很可能偏得很厉害。我自己一般会至少跑 10000 次并且在论文或者技术报告里注明是“蒙特卡洛经验门限”。4.4 我在复现时踩过的三个坑第一个坑复数噪声功率归一化。一开始我用randn(N,M)直接生成复噪声忘记除以 (\sqrt{2})结果噪声功率是 2 而不是 1。这样门限虽然还是按同样的统计量标定但信噪比定义全乱了曲线横轴整体偏移看起来像性能变好其实是假的。第二个坑门限方向写反。我把判决条件写成了T gamma判 (H_1)结果检测概率曲线几乎贴着 (P_{fa}) 走怎么调 SNR 都上不去。后来才反应过来(H_1) 下 (T_{PR}) 是变小而不是变大。这个错误很隐蔽因为你去看统计量的公式如果不明确写出每个假设下的趋势很容易搞反。第三个坑不同节点的信噪比设成完全一样忽略了近远效应。实际协作感知里距离主用户近的节点信噪比高远的节点信噪比低统一设成同一 SNR 会给出过于乐观的结果。严谨的仿真应该按路径损耗模型给每个节点分配不同的平均接收功率。PR检测器在那种情况下依然有效但检测概率曲线会比理想情况平缓不少。5. 工程落地中的扩展与个人建议5.1 从数据级融合到信息级融合集中式PR检测器的前提是融合中心能拿到所有节点的原始采样数据。这个前提在真实系统里很贵因为上传原始采样需要占用回传链路的大量带宽。工程上常用的折中方案是节点端先做降采样或压缩把观测数据压缩成若干个和统计量再传给融合中心。有人做过这样的尝试每个节点先计算自己接收信号与本地参考之间的相关值或者直接上传能量值融合中心把上传值拼成向量再构建协方差矩阵。这种做法的本质是把数据级融合降级成特征级融合性能会有损失但开销小得多。PR检测器在这种中间架构下依然可以使用只要注意协方差矩阵的维度是节点的数量而不是采样点的数量。5.2 多模态感知数据融合的一个方向最近总能看到“多模态感知数据融合”这个词。其实PR指数的核心思路——用协方差结构判断目标是否存在——也可以迁移到多模态感知里。假设你有温度、振动、射频等多个模态的传感器把各模态的特征拼成一个矩阵当某个目标事件发生时模态之间会出现统计相关性PR指数就能用来做事件检测或质量评估。不过跨模态有一个前提要小心不同模态的量纲和动态范围差异很大直接拼矩阵会把方差大的模态主导整个协方差结构。实际应用中要先把各模态数据做标准化或归一化再计算协方差。这个思路和我前面说“噪声功率会从统计量中抵消”是类似的逻辑完全一致。5.3 我的代码与参数使用建议最后给一点可以直接带走的建议。如果你想把蒙特卡洛跑得更快可以把最外层SNR循环改成parfor但要注意随机数流问题。Matlab的parfor默认会让每个worker从同一个全局随机流开始这会导致每个SNR点上使用的随机数完全一样破坏统计独立性。稳妥的做法是在循环内部用randn前显式设置不同的随机种子或者用RandStream给每个迭代分配独立的子流。离线门限标定和在线检测可以拆成两步。门限标定只需要做一次标定好之后存成gamma.mat在线检测时直接加载门限每次只需要计算一次prStat(Y)耗时几乎可以忽略。这个思路在硬件原型验证里也非常实用。我在复盘这个项目时最大的体会是PR检测器不是一个“性能碾压一切”的神器它真正的价值在于把“依赖绝对噪声功率”变成“依赖特征值相对结构”从而让频谱感知在噪声不确定度下依然可控。你只要把门限标定和理解每个假设下的统计量趋势这两件事做扎实复现这个检测器就非常顺手。仿真代码写完以后建议你自己改一改 (N)、(M)、(P_{fa})多画几张曲线对比一下很多直觉性的理解就会从这些对比里冒出来。