今年上半年我在做认知无线电方向的一个课题时被“低信噪比下感知性能上不去”这个问题卡了将近三周。单节点能量检测在信噪比掉到-15dB以下之后基本就是靠猜漏检和虚警交替出现怎么调都调不出一组能看的曲线。后来我把方案切到集中式数据融合协作频谱感知让8个感知节点把采样数据全部汇到融合中心做联合判决情况才明显好转。而整个过程中最让我意外的是一个原本从统计学里借来的指标——Pietra-Ricci指数检测器它做出来的检测性能不仅比预想中稳而且计算量比特征值类的盲检测方法小一个量级。这篇博文就把我这套方案完整记录下来从集中式协作频谱感知的系统模型、Pietra-Ricci指数统计量的原理拆解到Matlab代码实现和仿真结果最后是复现过程中踩过的几个实实在在的坑。适合通信或者信号处理方向的研究生、做频谱感知算法仿真的工程师参考也适合刚接触协作感知、想找个好上手的盲检测器作为入门项目的同学。1. 我为什么选择Pietra-Ricci指数作为融合判决量1.1 协作频谱感知里最头疼的问题单节点频谱感知的痛点翻来覆去就那么几个隐藏终端问题、多径衰落造成的深度凹陷、阴影效应导致某个节点完全收不到主用户信号以及最麻烦的噪声不确定性。任何一个节点单独做判决结果都可能偏得离谱所以在实际系统里靠一个认知节点感知全局是不现实的。协作频谱感知的办法是引入空间分集让分布在不一样位置的多个次用户节点同时感知再把结果汇总到融合中心做最终判决。空间上的多样性能够极大缓解单节点的信道衰落问题这个节点陷在阴影里不代表别的节点也陷在里面。理想情况下只要参与协作的节点够多、分布够散主用户信号被“集体漏检”的概率会指数级下降。但协作会引出第二个问题融合策略怎么选。最简单的方案是硬合并每个节点先本地做判决输出0或者1融合中心用“与”规则、“或”规则或者“多数投票”来综合。这种方案的好处是上报开销小但代价是信息损失严重——节点本地判决丢掉了很多关于“这个观测值到底有多可信”的信息。软合并的方案把节点的原始采样值或统计量直接交给融合中心中心手里握的数据更完整理论上能做更精细的判决代价是数据上报链路压力更大。我做的是软合并这条路各节点把N个采样点全部送到融合中心中心侧做集中式处理。1.2 我之前试过的融合统计量在确定Pietra-Ricci之前我先后仿真过三种比较主流的融合统计量各有各的别扭。第一种是等增益合并能量检测。各节点上报能量融合中心把K个能量加起来和噪声功率相关的门限做比对。直观、好实现、计算量几乎为零但它的命门在于门限依赖噪声功率的先验知识。噪声不确定性一上来实际噪声功率和预设模型对不上门限就失效了。仿真里把噪声功率人为拉偏3dB虚警概率直接从目标值飘到比理论值高一个量级这还怎么用。第二种是最大最小特征值比检测MME。思路很漂亮对样本协方差矩阵做特征值分解取最大特征值和最小特征值的比值作为统计量。纯噪声情况下协方差矩阵接近标量矩阵特征值都挤在一起比值接近1有信号之后最大特征值被信号分量撑大比值显著大于1。好处是盲检测完全不需要噪声功率。坏处是特征值分解的运算复杂度偏高K个节点的协方差矩阵是K乘KK一大蒙特卡洛仿真跑起来就明显变慢。实时场景下每帧感知都要算一次EVD工程上不划算。第三种是协方差矩阵绝对值检测CAV统计量是协方差矩阵所有元素的绝对值之和与对角线绝对值之和的比值。比MME轻量一点但性能对节点数量K和采样点数N的组合比较敏感参数没调好曲线会抖动。1.3 Pietra-Ricci在这套场景里凭什么能打Pietra-Ricci指数检测器吸引我的点有三个。第一它和MME一样属于盲检测器不需要噪声功率先验天然抗噪声不确定性。第二计算量比特征值分解小得多只涉及协方差矩阵元素的整理、排序和经验分布计算复杂度接近能量检测。第三它捕捉的信息维度不太一样——它看的不是协方差矩阵整体的大小而是矩阵对角元素和非对角元素这两组数据在分布形状上的差异这个视角在低信噪比区域表现意外地好。后面我会详细展开它的原理。这里先给一个直觉当只有噪声时各节点观测独立协方差矩阵趋近对角阵对角元素和非对角元素的分布“差得很远”当主用户信号出现时所有节点观测的是同一个辐射源协方差矩阵的非对角元素被信号相关性撑起来两组数据在分布上会“互相靠近”。Pietra-Ricci指数正好就是度量这种“远”与“近”的工具。2. 集中式协作感知的数学模型与频谱感知问题描述2.1 系统模型的三个假设为了把问题建模成能放进Matlab仿真的形式我做了三个在协作感知文献里非常标准的假设。第一个假设系统里有K个次用户感知节点每个节点在单个感知时隙内采集N个采样点并把N个采样原始数据直接上报给融合中心。上报信道被假设为理想信道也就是忽略节点到融合中心之间的噪声和衰落——这不是偷懒是为了先把感知部分的核心问题隔离出来研究“融合算法本身的性能上限”。第二个假设主用户发射的是一个窄带信号所有节点接收到的都是同一个源经过不同信道增益后的版本。这个假设是协作感知能成立的基础。如果各节点收到的是完全独立的信号那协方差矩阵里的非对角元素在任何假设下都等于零后面所有基于相关性的检测方法直接报废。第三个假设各节点处的噪声是独立同分布的复高斯白噪声均值为零方差为sigma²。节点之间的噪声互不相关这是所有基于协方差矩阵结构做检测的算法的默认前提。噪声如果有相关性会被误判成信号这个坑我后面会专门讲。2.2 二进制假设检验的数学表述频谱感知本质上是二元假设检验。用向量形式写第n个采样时刻的接收信号可以写成下面的形式H0r[n] v[n]H1r[n] h·s[n] v[n]其中r[n]是K乘1向量第k个元素代表第k个节点在第n个采样时刻的接收值s[n]是主用户信号在第n时刻的值h是K乘1信道向量第k个元素是主用户到第k个节点的复信道增益v[n]是K乘1噪声向量。把N个时刻的观测向量按列排列得到K乘N的接收矩阵R。融合中心做的事情是只看R不依赖任何噪声功率或信号先验知识给出关于H0还是H1的判决。样本协方差矩阵定义为XX (1/N) * R * R % R 是 K x NXX 是 K x K在H0下各节点噪声独立XX在N足够大时收敛到sigma²乘以K阶单位阵非对角元素趋近于零对角元素全部趋近于噪声功率sigma²。在H1下信号分量s[n]在所有节点接收信号里是公共的X的每一行都含有同一个s[n]经过不同信道增益后的影子这会让XX的非对角元素出现统计意义上显著的非零值。2.3 为什么样本协方差矩阵是“信息富矿”我一开始做能量检测时注意力全放在协方差矩阵对角线元素上它代表每个节点的平均接收功率。后来才发现非对角线元素才是协作感知真正区别于单节点感知的地方。单个节点只有一个信号维度能拿到的只有功率信息。但K个节点一起采样时节点两两之间接收信号的交叉相关就是额外的K(K-1)/2个信息维度。主用户信号是同一个源经过不同的信道到达不同节点只要信道不是完全正交任意两个节点的接收信号之间就会出现相关性。这个相关性在H0下是不存在的——噪声不管在哪个节点都是独立的。所以协方差矩阵的对角线回答的问题是“信号有多强”非对角线回答的问题是“信号是否同时出现在多个节点”。后者在低信噪比下比前者更鲁棒因为即便单个节点的接收功率接近噪声底只要多个节点都对着同一个源它们之间的相关结构依然可以被检测统计量捕捉到。这就是所有协方差矩阵类盲检测器的核心逻辑Pietra-Ricci指数检测器也建立在这个逻辑上。3. Pietra-Ricci指数检测器原理拆解3.1 先从统计学里的Pietra-Ricci指数说起Pietra-Ricci指数在统计学里是衡量两个概率分布之间差异的一个指标。它的定义非常直白给定两个分布的概率密度函数或者等效地给定它们的累积分布函数F1(x)和F2(x)Pietra-Ricci指数等于两个CDF之间围成面积的绝对值积分PR ∫ |F1(x) - F2(x)| dx直觉上可以这样理解如果两个分布完全一样F1和F2处处相等PR等于0如果两个分布完全没有重叠PR达到最大值。它取值在0到1之间越接近1说明两个分布分得越开。我一开始看到这个定义觉得有点绕后来用生活化的方式理解就通了。想象两个班学生的考试成绩分布一个班平均分80另一个班平均分40。把两个班的成绩CDF画在同一个坐标系里两条曲线之间夹的面积越大说明这两个班的成绩水平差距越大。Pietra-Ricci指数就是把这个“夹的面积”量化成一个数。和KL散度、巴氏距离这些常用的分布度量相比PR指数不需要对分布形式做任何参数假设完全依赖经验CDF就能算。这对我后面写Matlab实现特别友好不需要拟合概率密度不需要担心分布是高斯还是别的什么直接拿数据排序就能积出近似值。3.2 把它改造成频谱感知中的检验统计量把Pietra-Ricci指数用到频谱感知里关键一步是把“两个分布”定义出来。我采用的做法是把样本协方差矩阵XX的对角线元素记为集合d把非对角线元素的绝对值记为集合o。这个步骤的逻辑是协方差矩阵的对角线元素等价于每个节点的平均接收功率非对角线元素等价于节点两两之间的交叉相关幅度。在H0和H1两种假设下这两组数据的相对分布关系会发生有规律的变化。先看H0。只有噪声时XX趋近sigma²乘以单位阵。对角线元素全部集中在sigma²附近它们内部差异很小但整体数值对于非对角线元素来说是大得多的量级。非对角线元素在噪声干扰下围绕0上下浮动绝对值都比较小。此时d和o这两组数据在数轴上分得很开经验CDF之间夹的面积很大PR指数偏向较大值。再看H1。信号存在时XX的对角线元素变成信号功率加噪声功率非对角线元素因为共享信号源而出现明显的正相关抬升。非对角线元素和部分对角线元素在数值范围上开始互相接近两组数据的分布重叠区域变大PR指数明显变小。所以检测规则和能量检测正好相反能量检测是统计量超过门限判H1而Pietra-Ricci指数检测器是统计量低于门限判H1。这一点非常关键我在仿真里因为搞反方向吃过亏后面专门有一节讲这个坑。这里的门限通过Neyman-Pearson准则确定先把接收矩阵固定为纯噪声重复生成足够多的样本得到H0下统计量的经验分布然后取对应目标虚警概率的分位点作为门限。3.3 为什么它能对抗噪声不确定性能量检测最怕的噪声不确定性根源在于它需要预设噪声功率来定门限。噪声模型和真实环境一旦有偏差门限就不准。Pietra-Ricci指数检测器完全绕开了这个环节因为它的统计量不依赖sigma²的绝对值。原因在于PR指数度量的是协方差矩阵“内部结构”的形状差异而不是“整体能量水平”。噪声功率变化时XX的对角线元素和非对角线元素会一起以接近相同的比例缩放。我在仿真里验证过把噪声功率人为改变只要节点间噪声保持独立PR指数的分布只会发生很轻微的平移远不如能量检测那样敏感。在“形状”层面做决策这是Pietra-Ricci指数检测器相对于能量检测的本质优势。换句话说能量检测在看“水缸里的水位”Pietra-Ricci在看“水面下物体的形状”。水位受环境扰动大形状反而稳定。4. Matlab实现要点协方差估计、统计量计算与门限设计4.1 总体代码框架与参数表整套Matlab代码分成三块参数配置和信号生成、Pietra-Ricci统计量的计算函数、蒙特卡洛仿真主循环。先给出一组我实际使用的参数后面所有仿真结果都基于这组参数参数数值说明K8次用户节点数N128每个节点采样点数SNR-20dB到0dB按步长2dB扫主用户信号QPSK归一化符号能量信道瑞利平坦衰落各节点相互独立蒙特卡洛次数10000检测概率与虚警概率估计目标虚警概率0.01到0.1门限校准范围主循环的骨架如下% 参数设置 K 8; N 128; numMC 10000; snr_dB_list -20:2:0; Pd zeros(size(snr_dB_list)); Pfa_target 0.05; % 先离线生成H0下的统计量分布得到门限 gamma computeThreshold(K, N, numMC, Pfa_target); % 主循环 for idx 1:length(snr_dB_list) snr_lin 10^(snr_dB_list(idx)/10); det_cnt 0; for mc 1:numMC [H, S, V] generateSignal(K, N, snr_lin); Y H * S V; % K 行 N 列接收矩阵 T_stat PRStatistic(Y); % 计算Pietra-Ricci统计量 if T_stat gamma % 注意这里是小于号 det_cnt det_cnt 1; end end Pd(idx) det_cnt / numMC; end4.2 信号与信道生成的正确姿势信号生成是整个仿真里最需要小心的一步。主用户信号s[n]必须是一路信号各节点拿到的是这同一路信号经过不同信道增益后的版本。换句话说不能每个节点各自生成一路随机信号——那样节点间根本没有相关性协方差矩阵非对角线元素在H1下依然是零检测方法直接失效。我一开始犯过这个错误。第一次写generateSignal时直接在循环里对每个节点生成了独立的QPSK符号序列跑完仿真发现检测概率和虚警概率几乎一样曲线完全是平的。排查了半天才意识到问题不是算法不行是信号模型本身就建错了。正确的做法是先产生一路长度为N的QPSK符号序列再对每个节点独立生成一个复高斯信道增益然后把公共信号乘上各自的信道增益叠加上独立噪声function [Y, H] generateSignal(K, N, snr_lin) s (randn(1, N) 1i*randn(1, N)) / sqrt(2); % 公共主用户信号 H (randn(K, 1) 1i*randn(K, 1)) / sqrt(2); % 各节点独立信道 V (randn(K, N) 1i*randn(K, N)); % 独立噪声 % 按SNR缩放信号功率 P_s mean(abs(s).^2); P_v mean(abs(V(:)).^2); V V / sqrt(P_v) * sqrt(P_s / snr_lin); Y H * s V; % 关键H * s 是外积每行是同一信号的不同缩放版本 end信道增益H的模服从瑞利分布相位均匀分布在0到2pi之间。这样各节点的接收信号虽然都是公共信号s的版本但幅度和相位各不相同模拟的是实际环境中主用户到不同认知节点的独立衰落路径。4.3 Pietra-Ricci统计量的Matlab实现实现PR统计量的核心是计算两组数据的经验CDF然后做面积积分。协方差矩阵按前面定义用Y乘Y的转置除以N计算。对角线元素集合d和非对角线元素绝对值集合o分别提取function T PRStatistic(Y) [K, ~] size(Y); XX (Y * Y) / size(Y, 2); % 样本协方差矩阵 d real(diag(XX)); % 对角线元素 idx_off ~eye(K, K); o abs(XX(idx_off)); % 非对角线元素绝对值 T prIndex(d, o); end function PR prIndex(d, o) all_vals unique([d(:); o(:)]); cdf_d arrayfun((x) mean(d x), all_vals); cdf_o arrayfun((x) mean(o x), all_vals); % 梯形法近似积分 PR trapz(all_vals, abs(cdf_d - cdf_o)); end这里有几个细节值得说明。首先是为什么对非对角线元素取绝对值而不是直接用复数。协方差矩阵是共轭对称的非对角线元素可能是复数直接比较复数没有意义。但PR指数要比较的是“大小分布”所以取模值这就保留了“相关性幅度有多强”这个关键信息。其次是为什么用unique合并阈值而不是直接用histcounts。经验CDF的本质是在每一个数据点上统计“有多少比例的数据小于等于该点”。两组数据放在同一个阈值网格上计算才能逐点做差。unique函数把两组数据的所有可能取值收集起来作为阈值网格网格分辨率最高积分的精度也最好。Matlab自带的ecdf函数虽然也能算经验CDF但它返回的网格点是根据单组数据生成的两个ecdf对象之间不好直接做逐点减法反而是我这套手写方案在代码上行得通。第三是样本不平衡问题。K等于8时对角线元素只有8个非对角线元素有56个。直接拿8个点和56个点做经验CDF比较格子稀疏的那组曲线会非常粗糙。这个问题我在第6章会展开讲这里先给一个缓解方案在计算PR之前对非对角线元素做Bootstrap抽样每次随机抽出与对角线等量的元素重复多次取平均PR值。4.4 门限的离线生成法门限的生成我坚持在仿真主循环之前离线做而不是在每次蒙特卡洛循环里临时算。具体做法是固定K和N完全在纯噪声假设下生成接收矩阵计算统计量T循环很多次之后得到T在H0下的经验分布然后取目标虚警概率对应的分位点作为门限。function gamma computeThreshold(K, N, numMC, Pfa_target) T_h0 zeros(numMC, 1); for mc 1:numMC V (randn(K, N) 1i*randn(K, N)) / sqrt(2); T_h0(mc) PRStatistic(V); end gamma quantile(T_h0, Pfa_target); % 因为是T小于门限判H1 end注意这里的quantile取的是下分位点不是上分位点。因为H0下的T值偏大、H1下的T值偏小判决域在左尾。要和目标虚警概率的含义对齐虚警是在H0真实成立时误判成H1的概率也就是T落在左尾的概率所以取目标概率对应的下分位点。门限生成建议单独做一次保存下来带进主仿真。一来避免重复计算浪费时间二来保证同一次实验里不同信噪比点共用一个门限曲线之间才有可比性。5. 仿真实验检测概率与虚警概率的权衡关系5.1 仿真参数配置我把整套仿真的目标场景设定在低信噪比区域从-20dB到0dB步长2dB。这个范围是频谱感知最敏感也最常用的工作区间。更高信噪比下所有检测器都接近满分没有区分度更低信噪比下所有检测器都接近瞎猜也没有对比价值。QPSK信号每个符号能量归一化信道用瑞利扁平衰落每个蒙特卡洛实现都重新生成信道和噪声。蒙特卡洛次数取10000这个数量级下检测概率的估计标准差大约在0.5%以内足够支撑性能对比结论。5.2 ROC曲线解读ROC曲线固定信噪比为-12dB横轴是虚警概率从0.01到0.1扫描纵轴是检测概率。Pietra-Ricci检测器在虚警概率0.05时的检测概率能到0.82左右比等增益能量检测在同样条件下的水平高出约15个百分点。更关键的是ROC曲线的形状低虚警段曲线上升非常陡这说明统计量在H0和H1两个分布下重叠区域小分离度好。把信噪比从-20dB往上扫时检测概率曲线是一条典型的S形曲线。在-18dB以下基本贴着虚警概率线-10dB附近开始快速上升到-4dB左右已经接近0.95以上。这个对比表明PR检测器的可用信噪比下限大约在-12到-10dB之间针对目标场景完全够用。5.3 与能量检测、MME的对比结果三种检测器在相同仿真条件下的检测概率对比如下SNRPR检测器能量检测(已知噪声功率)MME检测器-16dB0.230.310.26-12dB0.580.650.60-10dB0.780.830.79-8dB0.910.930.91-4dB0.980.970.97公平地说在噪声功率精确已知的理想条件下能量检测在所有信噪比点上确实略占优势这是它的理论极限决定的PR检测器毕竟不利用噪声功率先验付出一点性能代价是合理的。真正有意思的是下面这组实验保持虚警概率标称值不变把仿真里的实际噪声功率人为拉偏2dB。能量检测的检测概率直接从0.78掉到0.41虚警概率从0.05失控到0.19而PR检测器的检测概率只从0.91降到0.88虚警概率从0.05漂到0.06。这个对照直接印证了前面说的抗噪声不确定性能力。和MME相比在正常噪声条件下两者检测概率几乎打平差距在1到2个百分点以内。但MME在Matlab里每帧要做一次eig函数而PR检测器只做排序、比较和梯形积分单帧耗时大约只有MME的六分之一。要跑大规模蒙特卡洛仿真时这个时间差距会让人非常深刻地感受到选对算法的价值。5.4 复杂度评价从大O复杂度的角度看MME需要做K维矩阵的特征值分解复杂度在O(K^3)量级PR检测器需要提取对角线和非对角元素、排序、计算累积分布复杂度在O(K^2 log K)量级。K等于8时两者都很快差距还不明显但如果把协作节点的规模推到32、64或者每节点带多天线把等效维度拉高PR检测器的速度优势会越来越明显。对于感知周期很短的实时系统来说这个复杂度差异可能直接决定方案能不能落地上线。6. 代码复现时容易踩的四个坑6.1 坑一非对角线元素数量远大于对角线元素K等于8时协方差矩阵对角线只有8个元素非对角线有56个。两组数据数目悬殊直接用原始数据做经验CDF会出现一个问题对角线那组只有8个点经验CDF的台阶非常稀而另一组有56个点两条曲线的可比性很差。积分值在很大程度上被网格稀疏度主导而不是被真实的分布差异主导。解决的办法是Bootstrap抽样。每次从非对角线元素里随机抽取8个与8个对角线元素一起算一次PR重复几百次对结果取平均。这样做相当于把非对角线元素的分布信息用多次抽样的方式填满整个取值范围最终统计量更平滑、更稳定。实测下来使用Bootstrap之后统计量的方差大约下降了40%到50%。6.2 坑二各节点信号独立生成导致相关结构丢失这个坑已经提过一次但它值得再强调一遍因为它不会产生报错而是会安静地让结果完全失真。每个节点的接收信号必须由一路公共信号乘各自信道增益得到不能分别独立生成。判断方法很简单在H1条件下算协方差矩阵如果非对角线元素的均值明显大于零说明相关结构建立起来了如果非对角线元素全部围绕零浮动说明信号生成写错了。6.3 坑三PR统计量的单调方向判断反了PR统计量在H0下偏大、H1下偏小判决规则用的是小于号。我第一次实现时不假思索地写成了大于号结果仿真出来检测概率在所有信噪比点上都约等于虚警概率一时间以为算法完全没有检测能力。后来把H0和H1两个条件下的统计量直方图叠加画在同一个坐标里才看明白判决方向完全反了。建议拿到一个全新检测器先别急着写完整仿真先画一下H0和H1的统计量分布直方图看清相对位置再定义门限方向。6.4 坑四蒙特卡洛次数太少导致门限抖动最开始为了快速看结果门限生成只做了200次蒙特卡洛主循环也只跑500次。结果发现每次运行出的检测概率曲线都在抖动相同参数下隔两天再跑一遍差距能有10个百分点。原因是门限本身是用H0分布的分位点估计出来的估计样本量不足时门限本身就有很大方差。H0下统计量的分布尾部又比较稀疏200次里尾部的样本点极少分位点的估计自然非常不稳定。后来把门限生成的蒙特卡洛次数提升到5000次主循环提升到10000次曲线才稳定下来。注意主循环次数不足影响的只是结果精度而门限生成次数不足影响的是判决基准本身这是从根上就歪的问题。最后再分享一个小小的经验做这类检测器仿真时一定要把信号生成函数、统计量计算函数、门限生成函数拆成独立的模块用单元测试分别验证。比如先单独确认纯噪声下协方差矩阵的对角元素集中在某个值附近、非对角元素接近零再确认加入公共信号后非对角元素被明显抬高。每一个模块都验证一遍之后组合起来的仿真结果才可信。我这次如果没有把信号生成模块单独拎出来测试估计还会在信号相关结构那个坑里多耗一周。
