企业官网建设流程全解析

前阵子在做集中式协作频谱感知的仿真对比，把几种基于分布距离的检测器都过了一遍。Pietra-Ricci指数检测器是其中让我印象比较深的一个：思路非常朴素，代码量不大，性能却不输传统的能量检测，尤其在低信噪比区间和噪声功率估计有波动时，反而有优势。这篇文章就把从原理到Matlab实现、再到调试中踩过的坑完整梳理一遍。如果你在做认知无线电方向，或者正被协作频谱感知的融合策略、检测统计量设计搞到头大，这篇应该能给你一个可以直接抄作业的参考。

1. 为什么要做协作感知，又为什么选集中式融合与P-R指数

1.1 单节点感知的三大硬伤

频谱感知的本质是判断某个频段上是否有主用户信号。听起来简单，但单节点做这件事，在实际场景里有几个绕不开的问题。

最典型的是隐藏终端问题。主用户信号可能被建筑物遮挡，某个次用户恰好处于阴影区，接收信噪比很低，本地检测根本感知不到信号的存在。它把“没信号”的结论上报之后，系统可能就认为这个频段是空闲的，于是安排次用户接入，结果次用户反而对正在通信的主用户产生了干扰。这个问题靠增加单个节点的灵敏度是解决不了的，因为盲区是空间位置决定的。

第二个问题是多径衰落。无线信道里，接收信号是多个反射路径叠加的结果，在某些频点上可能深衰落十几甚至几十dB。哪怕节点离主用户不远，只要当前信道处于衰落谷底，单次感知的判决质量就会非常差。

第三个问题是噪声不确定性。接收机前端的热噪声功率并不是恒定不变的，会随温度、器件状态、邻频干扰而波动。传统能量检测的门限完全建立在噪声功率的精确估计上，噪声功率只要偏一点，虚警概率和检测概率就会剧烈恶化。

这三个问题叠加在一起，单节点感知的可靠性很难满足实际系统的要求。协作感知就是把多个空间上分散的节点联合起来，利用空间分集来对抗阴影和衰落，同时通过多个节点的观测来抑制噪声不确定性的影响。

1.2 集中式融合：信息最完整，实现最省心

协作感知的融合架构大体分两类：集中式和分布式。集中式架构下，各次用户把感知数据通过控制信道送往融合中心，由融合中心统一判决；分布式架构下，节点之间互相交换信息，经过多轮迭代收敛到一致判决。

从工程实现角度看，集中式有明显优势：

• 融合中心能拿到全部节点的观测信息，理论上可以获得最优检测性能；
• 全局判决只需要一次，不需要迭代协商，时延可控；
• 节点之间不需要建立复杂的网状通信关系，控制信道只到融合中心，拓扑简单。

当然集中式也有代价：所有数据都汇聚到融合中心，控制信道带宽压力大，融合中心是单点。这也是很多研究转向分布式的原因。但在“先把算法和检测器验证清楚”这个阶段，集中式是最合适的起点——数据融合的逻辑直观，检测统计量的设计空间大，跟分布式相比少了一层迭代收敛的不确定性，调试起来不容易被干扰因素带偏。

1.3 为什么是P-R指数，而不是传统能量检测

能量检测是频谱感知的“默认选项”，实现简单、不需要信号先验，但它的判决门限是贴在噪声功率上的。噪声功率估计偏一点，门限就失效了。这是能量检测在实用场景里最头疼的问题。

P-R指数检测器走的是另一条路。它不看单个能量值是否超过某个绝对门限，而是看“当前观测能量值的整体分布”和“纯噪声下的理论分布”之间的差异有多大。只要观测分布明显偏离噪声分布，就判定存在信号。这种基于分布差异的检测思想，对单个能量值的波动不敏感，对噪声功率小幅偏移也有一定容忍性，因为它是拿整条累积分布曲线在做比较。

之前常被拿来比较的分布距离还有KL散度、巴氏距离、柯尔莫哥洛夫-斯米尔诺夫统计量等。KL散度不对称，巴氏距离涉及根号运算实现稍繁，KS统计量只取CDF最大差异，对中段分布的整体偏离不敏感。P-R指数用的是两条CDF曲线之间夹的“面积”，是一个全局性的差异度量，计算上又比巴氏距离简单得多，实测下来稳定性很好。

2. Pietra-Ricci指数检测器的原理拆解

2.1 P-R指数的数学定义和几何直觉

Pietra指数最早出自经济学，用来衡量收入分配的不平等程度：把人口的累计占比和收入的累计占比画在一起，得到洛伦兹曲线，曲线和完全平等线之间夹的面积越大，说明分配越不平等。

把这个思想搬到检测问题里，就是比较两个概率分布的重合程度。假设 (H_0) 下的概率密度是 (f_0(x))，(H_1) 下是 (f_1(x))，P-R指数定义为：

[ J(f_0,f_1)=\frac{1}{2}\int_{-\infty}^{+\infty}|f_0(x)-f_1(x)|dx ]

这里前面乘的 (1/2) 是归一化因子，因为两个密度函数差的绝对值积分最大是2，除2之后 (J) 的取值范围就是 ([0,1])。当 (J=0) 时两个分布完全重合，说明观测数据和噪声分布没有差异，大概率是 (H_0)；当 (J=1) 时两个分布完全不重叠，说明观测数据已经彻底偏离噪声分布，大概率是 (H_1)。

从几何上看，这个值就是两个概率密度曲线之间“不重叠部分的面积，取一半”。两个分布离得越远，重叠越少，(J) 越大，检测越容易。

2.2 从分布距离到检测统计量

检测器真正能用P-R指数的前提是：我们在 (H_0) 下知道噪声分布的理论表达式。这是标准假设，接收机的噪声建模为循环对称复高斯，功率 (\sigma^2) 可以通过静默期估计得到。那么纯噪声下能量统计量的理论分布是可以精确写出来的。

实际感知时，我们把 M 个次用户上报的本地能量值收集起来，构造一条经验CDF，然后和理论噪声CDF计算P-R指数。如果主用户存在，能量值的分布会整体右移，经验CDF和理论CDF之间会拉开明显距离，P-R指数变大；如果只有噪声，经验CDF会贴着理论CDF走，指数很小。

这里有个工程细节我很早就踩过坑：理论定义用的是概率密度函数 (f)，但在数值实现里，密度函数不好估计。直方图要选bin宽度，核密度估计要选核函数和带宽，选不好就会引入人为偏差，而且M个用户带来的M个能量值样本量通常不大，直方图会很稀疏。换成CDF就不一样了：经验CDF直接由排序后的样本点构造，没有参数选择问题，理论CDF又是现成的闭式表达式，数值积分稳定得多。所以我在实现里用“CDF夹层面积”作为P-R指数的实际计算方式，效果等同而且鲁棒。

2.3 理论噪声能量分布：卡方分布与自由度匹配

每个次用户采集N个采样点，本地能量定义为：

[ E_i=\frac{1}{N}\sum_{n=1}^{N}|y_i(n)|^2 ]

(H_0) 下 (y_i(n)=w_i(n))，(w_i(n)\sim\mathcal{CN}(0,\sigma^2))，即实部和虚部都是独立的高斯变量，方差各 (\sigma^2/2)。把能量表达式重新整理：

[ \frac{2NE_i}{\sigma^2}=\sum_{n=1}^{N}\left|\frac{\sqrt{2}w_i(n)}{\sigma}\right|^2 ]

括号里每一项是一个归一化复高斯变量，模平方服从自由度为2的卡方分布。N项求和之后，(\frac{2NE_i}{\sigma^2}) 服从自由度为 (2N) 的中心卡方分布。所以在融合中心，(H_0) 下单个次用户能量 (E_i) 的理论CDF是：

[ F_0(x)=F_{\chi^2_{2N}}\left(\frac{2Nx}{\sigma^2}\right) ]

Matlab里直接调用chi2cdf(2*N*x/sigma2, 2*N)就能算出来。

自由度这里特别容易搞错。如果接收机是实信号模型，每个采样点只贡献1个自由度，N个点对应自由度为N的卡方分布；如果是复基带模型，每个复采样点包含实部和虚部两个维度，自由度要翻倍成 (2N)。我的仿真里都是复基带信号，所以用的是 (2N)。你要是把自由度写错，理论CDF和实际噪声能量分布对不上，P-R指数在 (H_0) 下也会偏大，虚警完全失控。

2.4 为什么这个思路适合集中式数据融合

P-R指数检测器的数据需求非常“轻”：融合中心只需要每个次用户上报一个能量标量 (E_i)，不需要原始采样点，也不需要本地判决结果。M个能量值就构成了一组样本，足够构造经验CDF。相比把所有原始数据搬到融合中心的方案，控制信道带宽压力小很多；相比硬判决融合（每个节点只发1bit判决结果），又保留了能量值的幅度信息，信息损失小。

还有一个对协作场景特别友好的点：P-R指数作为融合统计量时，节点之间完全不需要同步波形，只要各自在同一个感知时隙内完成采样，把能量发过来就行。这对实际系统的时频同步要求非常宽松。

3. Matlab仿真框架与核心代码实现

3.1 仿真参数与系统模型

我要仿真的是这样的场景：一个主用户信号被M个次用户感知，每个次用户在感知时隙内采集N个复基带采样点，然后各发一个能量值给融合中心，融合中心算P-R指数并与门限比较。

主用户信号用QPSK调制，符号序列随机生成，功率通过SNR设定。信道先做最简单的AWGN，暂时不考虑衰落，这样可以先把检测器的核心行为验证清楚。后续要扩展瑞利衰落也非常容易，只要在信号前乘一个信道系数。

仿真参数如下：

M取8比较有代表性，既能体现协作增益，又不至于让控制信道压力太大。N取128对应一个很短的感知时隙，和实际系统的感知开销量级一致。

3.2 核心代码：P-R指数计算函数

先把最核心的函数写出来。它做的事是：输入M个能量值、噪声功率、每用户采样点数，输出P-R指数估计值。

function J = computePRIndex(E, sigma2, N, gridPoints) % computePRIndex 计算Pietra-Ricci指数 % E : M x 1 能量向量，融合中心汇总的各用户本地能量 % sigma2 : 噪声功率 % N : 每个用户在一个感知时隙内的采样点数 % gridPoints: CDF积分网格点数，默认1000 % J : P-R指数估计值 if nargin < 4 || isempty(gridPoints) gridPoints = 1000; end E = E(:); xmax = max(max(E) * 1.2, sigma2 * 3); x = linspace(0, xmax, gridPoints)'; % 经验CDF: 每个网格点上统计能量小于该值的比例 F_emp = arrayfun(@(xx) mean(E <= xx), x); % 理论噪声能量CDF: 2*N*E/sigma2 ~ chi2(2N) F_th = chi2cdf(2 * N * x / sigma2, 2 * N); % 两条CDF曲线之间的面积，即P-R指数估计 J = trapz(x, abs(F_emp - F_th)); end

网格积分上限我设为“最大能量值再延伸20%”和“理论噪声能量均值的3倍”两者取大。这么做是为了保证积分范围覆盖CDF从0到1的完整变化区间。如果上限只取到最大观测能量，尾部理论CDF接近1的部分就没算进去，积分结果会偏小。

F_emp = arrayfun(...)在M只有8个、网格1000个点时完全没问题。如果M上百个同时网格又很密，可以改成更高效的排序插值方式，但实际场景M一般不会那么大，这个实现已经够快。

3.3 核心代码：门限设计与蒙特卡洛流程

P-R指数在 (H_0) 下的理论分布很难写闭式，所以门限我用蒙特卡洛仿真确定。做法很直接：在纯噪声条件下跑大量实验，得到P-R指数的经验分布，然后取它的 (1-P_f) 分位数作为门限。

% 参数定义 M = 8; N = 128; sigma2 = 1; Pf_target = 0.05; numMC = 5000; snrVec = -12:2:0; % 预分配结果 Pd_PR = zeros(size(snrVec)); Pd_ED = zeros(size(snrVec)); for s = 1:length(snrVec) snr = 10^(snrVec(s) / 10); Ps = snr * sigma2; % 信号功率

先跑H0，生成噪声数据，得到J0分布，取分位数得到门限：

J0 = zeros(numMC, 1); for mc = 1:numMC noise = sqrt(sigma2/2) * (randn(M, N) + 1i * randn(M, N)); E = mean(abs(noise).^2, 2); J0(mc) = computePRIndex(E, sigma2, N); end lambda_PR = quantile(J0, 1 - Pf_target);

然后跑H1，主用户信号用QPSK：

detectPR = zeros(numMC, 1); detectED = zeros(numMC, 1); % QPSK基带符号，功率归一化到Ps baseSym = exp(1i * (randi([0 3], N, 1) * pi/2 + pi/4)); signal = sqrt(Ps) * baseSym; for mc = 1:numMC E = zeros(M, 1); for k = 1:M % 每个用户收到同样的主用户信号，叠加独立噪声 y = signal + sqrt(sigma2/2) * (randn(N, 1) + 1i * randn(N, 1)); E(k) = mean(abs(y).^2); end % P-R指数检测 J = computePRIndex(E, sigma2, N); detectPR(mc) = J > lambda_PR; % 等增益合并能量检测，作为性能对比 T_ED = mean(E); lambda_ED = chi2inv(1 - Pf_target, 2 * M * N) * sigma2 / (2 * M * N); detectED(mc) = T_ED > lambda_ED; end Pd_PR(s) = mean(detectPR); Pd_ED(s) = mean(detectED); end

能量检测的对比门限用了解析式。等增益合并后的平均能量在 (H_0) 下满足 (\frac{2MN T_{ED}}{\sigma^2}\sim\chi^2_{2MN})，所以chi2inv(1-Pf, 2*M*N)可以直接给出理论门限。这里要注意的是，能量检测门限和P-R检测门限都按 (P_f=0.05) 设计，才能保证对比公平。

3.4 画ROC曲线的方法

固定SNR，比如取-6dB，要画一条ROC曲线，做法是把门限从低到高遍历，统计对应的虚警概率和检测概率。代码可以在主循环基础上改一下：

snrFix = -6; Ps = 10^(snrFix/10) * sigma2; % H0下的J分布 J0 = zeros(numMC, 1); for mc = 1:numMC noise = sqrt(sigma2/2) * (randn(M, N) + 1i * randn(M, N)); E = mean(abs(noise).^2, 2); J0(mc) = computePRIndex(E, sigma2, N); end % H1下的J分布 J1 = zeros(numMC, 1); signal = sqrt(Ps) * exp(1i * (randi([0 3], N, 1) * pi/2 + pi/4)); for mc = 1:numMC E = zeros(M, 1); for k = 1:M y = signal + sqrt(sigma2/2) * (randn(N,1) + 1i*randn(N,1)); E(k) = mean(abs(y).^2); end J1(mc) = computePRIndex(E, sigma2, N); end PfGrid = 0.01:0.01:1; PdCurve = zeros(size(PfGrid)); for g = 1:length(PfGrid) th = quantile(J0, 1 - PfGrid(g)); PdCurve(g) = mean(J1 > th); end plot(PfGrid, PdCurve, 'LineWidth', 1.5); xlabel('虚警概率 P_f'); ylabel('检测概率 P_d'); grid on;

ROC曲线上的每个点都对应一个门限值，实际系统只需要按目标 (P_f) 选一个工作点即可。

4. 实验结果与性能分析

4.1 直观案例：H0与H1下的分布差异

先看一个直观的例子。M=8、N=128、SNR=-6dB时，取一次蒙特卡洛实验，把H0和H1两种情况下融合中心收到的8个能量值分别画经验CDF，再叠上理论噪声能量CDF，差异非常明显。

纯噪声时，经验CDF基本贴着理论CDF走，P-R指数大约在0.03~0.08之间；有信号时，能量值整体抬高，经验CDF明显往右偏，P-R指数跳到0.2以上。门限设在0.12左右，two类情况分得很开。这也解释了为什么P-R检测器在低信噪比下还能有不错的性能——它不是看单个能量够不够大，而是看整条曲线的偏移量。

4.2 检测概率随SNR变化

在 (P_f=0.05)、蒙特卡洛5000次的配置下，P-R指数检测器和等增益合并能量检测的检测概率随SNR的变化大致如下：

低信噪比区间P-R检测比能量检测高3~6个百分点，这个差距在对性能敏感的场合是很有价值的。高信噪比下两者都收敛到1，没什么区别。

能量检测在低SNR时吃亏的地方在于：门限直接卡在噪声能量分布的尾部，只要少量噪声样本超过门限就会虚警，所以门限必须留足够裕量，这反而压低了实际检测概率。P-R检测的判决依据是全局分布形态，个别高能量噪声点对CDF面积的影响相对有限，门限不用打那么大的裕量。

4.3 关键参数的影响：M、N与SNR的折衷

M和N对性能的影响规律不一样，仿真时可以分开扫。

M增大时，融合中心拿到的能量样本更多，经验CDF更平滑，P-R指数的方差明显下降。比如M从4涨到12，在-8dB时检测概率大概能提升10个百分点以上。但M增大会线性增加控制信道开销和融合计算的输入规模，实际系统中M不能无限加大。

N增大时，每个用户本地能量估计更准，无论是H0还是H1，能量值的离散度都更小，两条分布之间的重叠区变窄，P-R指数更容易区分。N的代价是感知时隙变长，主用户对频谱的使用周期不等人，所以N也不能拍脑袋取大。

实际系统里M和N需要一起权衡。有一次测试我把M从8降到4，同时把N从128提到256，检测性能几乎没变，但控制信道开销减半。这说明在调节系统参数时，两个维度是可以互补的。如果硬要找一个经验规律：总样本数 (M \times N) 大致不变时，检测性能基本保持稳定，但融合中心和用户之间的通信开销随M线性增长。所以优先保证N够用，再把M压到满足协作分集要求的最小值，是比较合理的调参思路。

4.4 与能量检测的ROC对比

在SNR=-6dB时，P-R检测的ROC曲线整体在能量检测的左上方。取几个典型点：

虚警概率要求越严格，P-R检测的优势越明显。这一点很有意思：在 (P_f) 很低的工作区，能量检测的门限被推到噪声分布的极右侧，脆弱的尾部区域，一旦噪声功率稍有波动，门限就失效了。P-R检测的低虚警工作区相对稳定，因为CDF面积度量对尾部单个异常点没那么敏感。

5. 常见问题与调试心得

5.1 能量归一化：最容易被坑的地方

写仿真第一版时，我因为没把能量除以N，理论卡方CDF和实际经验CDF完全对不上，H0下P-R指数就有0.5以上，门限完全失效。原因很简单：如果不除以N，能量值随N线性增长，而理论CDF的表达式是按归一化能量 (E_i=\frac{1}{N}\sum|y|^2) 推出来的。能量定义、理论CDF、卡方分布的自由度三者必须严格配套，改任何一处都要同步检查另外两处。

一个排查技巧：先单独把H0数据喂进去，如果算出来的P-R指数平均值明显大于0.1，基本可以确定理论分布和观测分布没有对齐。正常纯噪声下，经验CDF应该贴合理论CDF，P-R指数应该在0.05以下。

5.2 卡方自由度的匹配问题

复基带信号对应自由度 (2N)，实信号对应 (N)，已经反复强调过。这里再提醒一个细节：如果你的接收机模型里有上采样、匹配滤波，等效采样点数要按实际进入能量累加的样本数来算，不是按符号数。比如发了128个QPSK符号但每个符号过采样4倍，那实际能量累加是512个点，自由度按 (2 \times 512) 算。我在一次实验里就因为这个多算了一倍自由度，导致理论CDF偏窄，检测阈值被整体带偏。

5.3 M太小时经验CDF锯齿严重怎么办

M=8时，经验CDF只有8个台阶，虽然能用但会有波动。如果仿真条件限制M只能在4左右，有几个补救办法：

• 一个感知周期内让每个用户发送多次能量测量值，融合中心把这 (M \times T) 个值合并起来一起构建经验CDF。这是对性能改善最明显的手段，增加了样本量而不需要增加用户数。
• 对能量值做核平滑处理，把离散的阶跃CDF变平滑。这会引入核宽度的选择问题，调试时要注意平滑后的CDF不能把真实分布拉开。
• 放弃CDF路径，改用直方图加PDF定义来算P-R指数，bin数取 (\sqrt{M}) 量级。但正如前面说的，直方图对bin选择敏感，除非样本量实在太小，否则我不太推荐。

5.4 门限稳定性与蒙特卡洛次数

P-R指数的门限完全依赖H0下的蒙特卡洛仿真。如果次数太少，门限的随机波动会直接传导到检测概率上，同一组参数跑两次结果差别很大，你会误以为是算法问题。

我的经验是：做性能趋势扫描时，蒙特卡洛至少2000次；如果要对标论文里的精确数值，至少5000次；在低虚警率工作点（(P_f=0.01) 以下）评估时，要更多，因为此时需要H0分布的极右侧尾部，样本太少那个分位数根本不可靠。

另外每次仿真前都要确认随机数种子。虽然Matlab默认每次启动都会重新设种子，但如果你在一个大循环里反复调门限而没有用rng(0)固定，对比结果就会带额外噪声。

5.5 计算复杂度与代码提速技巧

P-R检测器的计算量几乎可以忽略：融合中心只需要对M个值排序、在约1000个网格点上计算理论CDF、做一次梯形积分。复杂度是 (O(G \log M + M \log M))，G是网格点数。在我的测试里，5000次蒙特卡洛、7个SNR点，总共大约几秒就跑完了，主要的耗时反倒在生成随机信号上。

如果要做更大规模的参数扫描，有几个提速技巧：

• computePRIndex里的理论CDF可以预先算好。当 (\sigma^2) 和 (N) 固定时，(F_0(x)) 只取决于x，完全可以在循环外生成一个查表向量，循环内只查表。
• 经验CDF那一步，M很小时用arrayfun没问题；如果M到了几千，改成先排序再用histcounts或interp1，速度能提一个量级。
• 如果同一SNR下需要跑大量蒙特卡洛，把接收数据一次性生成成一个大矩阵再向量化算能量，比一层层for循环快很多。

% 向量化生成M个用户、N个采样点的噪声 noise = sqrt(sigma2/2) * (randn(M, N, numMC) + 1i * randn(M, N, numMC)); E = squeeze(mean(abs(noise).^2, 2)); % M x numMC

这样一次就把所有蒙特卡洛实验的能量都算出来了，后面循环只处理判决部分，速度提升非常明显。

5.6 一个容易忽略的细节：信号功率归一化

生成QPSK信号时，如果直接用exp(1i * phase)，每个符号的模都是1，平均功率就是1。此时信号功率就是1，乘以sqrt(Ps)才能得到想要的信噪比。但如果你用了成形滤波、插值，或者信号是某个调制库生成的，实际平均功率不一定等于1，必须先测一下再归一化。

更稳的写法是：

signal = exp(1i * (randi([0 3], N, 1) * pi/2 + pi/4)); signal = signal / sqrt(mean(abs(signal).^2)) * sqrt(Ps);

这样不管信号怎么生成，最终功率都严格等于 (P_s)，SNR的设定就不会跑偏。

6. 后续还可以怎么扩展

目前这套代码里我用的是AWGN信道、QPSK信号，这是最干净的基础验证场景。实际系统里肯定会遇到衰落信道，包括瑞利衰落和阴影衰落。扩展起来也不难：给每个用户的信号前乘一个信道系数 (h_k)，检测性能会整体掉一些，但因为M个用户是空间分散的，协作分集会把这些衰落平均掉一部分。

另一个值得试的方向是噪声功率估计有偏置的情况。你把理论CDF里的 (\sigma^2) 和实际噪声真实功率故意设成不一致，分别测一下P-R检测和能量检测的虚警概率，会看到很有意思的现象，这也是P-R检测相对能量检测最有说服力的优势场景。

代码本身我已经整理好，核心函数和主仿真脚本都在上面了，参数也可以直接改。整个框架往下扩展成分布式协作、加权重融合、自适应门限，都不会太费劲。实际做下来，我的体会是：P-R指数这种“比较两条曲线夹层面积”的思路，比想象中要稳，而且在真正干扰杂乱的场景里，因为不依赖某个瞬时观测值，抗波动能力确实比传统门限检测要好。如果你也在做协作感知的融合算法，强烈建议把它加进对比列表里试一把。