基于Cholesky矩阵分解的协作频谱检测方法

摘要

本发明公开了一种基于Cholesky矩阵分解的协作频谱检测方法，主要解决现有的能量检测算法受噪声不确定影响大，需要先验信息和检测门限难以确定等问题。其实现步骤是：(1)各个检测用户根据所要检测的频段，采集该频段的数据，并上传到处理中心；(2)处理中心根据上传的数据，构造归一化协方差矩阵，并对其进行Cholesky分解；(3)利用该分解结果计算检测统计量，分析其概率分布；(4)根据检测统计量的概率分布计算目标虚警概率下的判决门限；(5)处理中心对检测统计量与检测门限进行比较，判决主用户信号是否存在。本发明具有抗噪声不确定度能力强，检测门限精准，检测性能高的优点，可用于无线通信。

说明书

技术领域

本发明属于无线通信技术领域，涉及一种频谱检测方法，可用于感知网络和认知无 线电系统中的频谱检测。

背景技术

随着感知网络的发展及其在日常生活中的普及，面向感知网络的关键技术的研究 得到了广泛关注。其中，频谱检测技术是感知网络中的重要关键技术，通过检测频段 的空闲与否来决定是否利用当前频段。同时，随着无线通信和移动通信的发展，人们 对通信技术的要求越来越高，所需的服务越来越多，这使得有限的频谱资源变得日益 稀缺。在目前的静态频谱分配框架下，很多频谱资源被分配给一些特定的服务，使得 频谱的利用率低下。为了改善频谱利用率低下的现状，J.Mitola等人提出了认知无线电 的概念，其主要思想是在已授权的频段内寻找空闲频段，在不影响授权用户正常通信 的前提下，允许检测用户能够检测并接入到当前空闲的频段，从而大幅提高频谱利用 率。为了能接入到空闲的频段，检测用户必须准确地检测其周围的频谱占用情况，因 此频谱检测技术在认知无线电中具有关键的作用。该项技术包括频谱检测方法和协作 频谱检测方法。

现有的频谱检测方法，主要有三种：

1)能量检测。检测用户通过计算接收到信号的能量，根据信号的能量大小来确定 主用户是否存在。该方法实现简单，容易确定检测门限。然而，在低信噪比情况下， 由于受到深度衰落和多径衰落等因素的干扰，该方法不能有效地正常工作。而且，该 方法受噪声不确定度的干扰，在实际应用中受限。

2)基于循环平稳的检测。检测用户利用主用户信号在循环频率处所表现的峰值特 性与噪声在循环频率处没有峰值特性来检测主用户是否存在。该方法抗噪性能好但需 要主用户的先验信息且复杂度高。在认知无线电中其系统效率低，在实际应用中受限。

3)基于特征值分解的检测。检测用户利用主用户信号的相关性，通过对协方差矩 阵进行特征值分解，来构造检测统计量。该方法能够抵抗噪声不确定度问题，且性能 优于能量检测。然而，该方法只能采用无限采样点以确定近似的检测门限，其检测性 能随之降低。

协作频谱检测是各检测用户通过协作来确定主用户存在与否。现有的协作频谱检 测方法，包括：

1)基于能量检测的协作检测。各检测用户利用能量检测方法，通过协作确定主用 户存在与否。该方法虽然实现简单，但是易受噪声不确定的干扰，而且，在低信噪比 情况下，容易受深度衰落和多径衰落等因素的干扰，不能有效地正常工作。

2)基于循环平稳的协作检测。各检测用户利用基于循环平稳的检测方法，通过协 作确定主用户存在与否。该方法抗造性能好，但复杂度高，且需要主用户信号的先验 知识，无法做到盲检测，在实际中受限。

3)基于特征值分解的协作检测。各检测用户利用基于特征值分解的检测方法，通 过协作确定主用户存在与否。该方法能够抵抗噪声不确定度的影响，检测性能好，但 难以精确地确定检测门限，给实际应用造成了局限。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于Cholesky矩阵分解的 协作频谱检测方法，以提高抗噪声不确定度能力，降低检测复杂度，精确确定判决门 限，提高对主用户信号的检测性能。

为实现上述目的，本发明的技术方法包括如下步骤：

(1)将占用当前频段的用户信号定义为主用户，将通过检测当前频段上主用户是 否存在，以试图占用该频段的用户信号定义为检测用户，将通过融合和分析各个检测 用户采集的数据，以确定当前频段主用户信号是否存在的设备定义为处理中心；

(2)各个检测用户根据所要观察的频段，采集该频段的数据得到x_i(n)，其中 n＝1,...,N；i＝1,...,M，N为各个检测用户的采样点数，M为检测用户数，各个检 测用户将采集到的数据x_i(n)上传到处理中心；

(3)处理中心根据各检测用户上传的数据x_i(n)，构建检测统计量T_ξ：

(3.1)处理中心根据各个检测用户上传的数据x_i(n)，构建数据矩阵X和协方差 矩阵R_x，其中数据矩阵X为：

$X=\left(\begin{array}{cccc}{x}_1\left(1\right)& x_1\left(2\right)& ...& x_1\left(N\right)\\ x_2\left(1\right)& x_2\left(2\right)& ...& x_2\left(N\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ x_M\left(1\right)& x_M\left(2\right)& ...& x_M\left(N\right)\end{array}\right),$

协方差矩阵为：

$R_x=\frac{1}{N}{\mathrm{XX}}^H,$

其中(·)^H为Heimitian转置；

(3.2)处理中心根据协方差矩阵R_x，计算归一化协方差矩阵R'_x：

$R_x^{\prime }=\frac{N}{{\sigma }_w^2}{R}_x,$

其中，为处理中心设置的噪声方差，N为各个检测用户的采样点数；

(3.3)处理中心归一化协方差矩阵R'_x进行Cholesky分解，得到分解后的上三角矩 阵，即：

R'_x＝L^TL，

其中，L为上三角矩阵，其表示为：

$L=\left(\begin{array}{cccc}{l}_{11}& l_{12}& ...& l_{1M}\\ 0& l_{22}& ...& l_{2M}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 0& 0& ...& l_{\mathrm{MM}}\end{array}\right),$

其中，l_ij为上三角矩阵L的第i行第j列元素，i＝1,...,M,j＝1,...,M；

(3.4)处理中心根据分解后得到的上三角矩阵L，构建检测统计量T_ξ：

$T_{\xi }=\frac{{\zeta }_1}{{\zeta }_M},$

其中ζ₁，ζ_M分别为上三角矩阵L的最大特征值和最小特征值；

(4)处理中心根据检测统计量T_ξ，计算检测门限γ_ξ：

${\gamma }_{\xi }=\frac{u_{{\zeta }_M}{u}_{{\zeta }_1}}{u_{{\zeta }_M}^2-{\Phi }^{-1}(1-P_{\mathrm{fa}}){\sigma }_{{\zeta }_M}^2}+\frac{{\Phi }^{-1}(1-P_{\mathrm{fa}})\sqrt{{\sigma }_{{\zeta }_1}^2{u}_{{\zeta }_M}^2+{\sigma }_{{\zeta }_M}^2{u}_{{\zeta }_1}^2-{\Phi }^{-1}(1-P_{\mathrm{fa}}){\sigma }_{{\zeta }_1}^2{\sigma }_{{\zeta }_M}^2}}{u_{{\zeta }_M}^2-{\Phi }^{-1}(1-P_{\mathrm{fa}}){\sigma }_{{\zeta }_M}^2},$

其中，和分别为步骤(3.4)中的最大特征值ζ₁和最小特征值ζ_M的均值，和分别为最大特征值ζ₁和最小特征值ζ_M的方差，P_fa为虚警概率，取值范围为 (0,1)，Φ^-1(·)为标准正态分布的累积量分布函数Φ(·)的逆函数；

(5)将步骤(3.4)中得到的检测统计量T_ξ与步骤(4)中得到的检测门限γ_ξ进行 比较，当T_ξ≥γ_ξ时，判决为主用户存在，即当前频段频谱已被某用户占用，否则，判 决为主用户不存在，即当前频段频谱为空闲状态，允许检测用户利用。

本发明具有以下优点：

1、本发明利用主用户信号的相关性进行检测，检测性能优于基于最大最小特征值 之比的协作检测方案以及基于最大特征值和矩阵迹之比的协作检测方案。

2、本发明是一种全盲检测方法，不需要任何有关主用户，信道和噪声的先验信息。

3、本发明基于有限随机矩阵的Cholesky分解，根据随机矩阵理论，得到检测门限 的闭式表达式，能在任意的采样点数下获得与目的虚警概率对应的精确的检测门限。

4、本发明能快速确定检测门限，降低了频谱检测复杂度，可以在实际中广泛应用。

附图说明

图1是本发明的实现流程图；

图2是本发明在主用户信号不存在情况下，针对理论和仿真，得到的归一化协方 差矩阵的最大特征值的概率分布的累积量分布的曲线对比图；

图3是本发明中在主用户信号不存在情况下，针对理论和仿真，得到的归一化协 方差矩阵的最小特征值的概率分布的累积量分布曲线对比图；

图4是本发明中在主用户信号不存在情况下，针对理论和仿真，得到的归一化协 方差矩阵的最大特征值和最小特征值的比值的概率分布的累积量分布曲线对比图；

图5是用本发明方法，现有的基于最大最小特征值之比的协作检测方法以及现有 的基于最大特征值和矩阵迹之比的协作检测方法，得到的“性噪比——检测概率”对 比图；

图6是用本发明方法，现有的基于最大最小特征值之比的协作检测方法以及现有 的基于最大特征值和矩阵迹之比的协作检测方法，得到的ROC曲线对比图。

具体实施方式

参照图1，本发明的实现步骤如下：

步骤1，各个检测用户采集数据，并上报给处理中心。

1.1)各个检测用户根据所要检测信号的频段，用相应的滤波器滤出该频段的信号；

1.2)在满足采样定理的前提条件下，对该频段的数据进行采集，得到采集到的数 据x_i(n)，其中n＝1,...,N；i＝1,...,M，其中N为各个检测用户的采样点数，M为 检测用户数；

1.3)将采集的数据x_i(n)上报给处理中心。

步骤2，处理中心根据上报的数据x_i(n)，得到归一化协方差矩阵R'_x。

2.1)处理中心根据各检测用户上报的数据，构建M行N列数据矩阵X：

$X=\left(\begin{array}{cccc}{x}_1\left(1\right)& x_1\left(2\right)& ...& x_1\left(N\right)\\ x_2\left(1\right)& x_2\left(2\right)& ...& x_2\left(N\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ x_M\left(1\right)& x_M\left(2\right)& ...& x_M\left(N\right)\end{array}\right),$

其中，N为各个检测用户的采样点数，M为检测用户数；

2.2)处理中心根据构建的数据矩阵X，计算协方差矩阵R_x：

$R_x=\frac{1}{N}{\mathrm{XX}}^H,$

其中(·)^H为Heimitian转置；

2.3)处理中心根据协方差矩阵R_x，计算归一化协方差矩阵R'_x：

$R_x^{\prime }=\frac{N}{{\sigma }_w^2}{R}_x,$

其中，为处理中心设置的噪声方差，取值为(0,+∞)。

步骤3，处理中心根据归一化协方差矩阵R'_x，按如下公式进行Cholesky分解：

R'_x＝L^TL，

其中，L为上三角矩阵，其表示为：

$L=\left(\begin{array}{cccc}{l}_{11}& l_{12}& ...& l_{1M}\\ 0& l_{22}& ...& l_{2M}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 0& 0& ...& l_{\mathrm{MM}}\end{array}\right),$

其中，l_ij为上三角矩阵L的第i行j列元素，i＝1,...,M,j＝1,...,M。

步骤4，处理中心根据上三角矩阵L计算检测统计量T_ξ。

处理中心根据概率论知识，计算得到上三角矩阵L的M个特征值为ζ₁,ζ₂,...,ζ_M， 其中ζ₁≥ζ₂≥...≥ζ_M，根据最大特征值和最小特征值计算检测统计量T_ξ，其表示为：

$T_{\xi }=\frac{{\zeta }_1}{{\zeta }_M},$

其中ζ₁，ζ_M分别为上三角矩阵L的最大特征值和最小特征值。

步骤5，处理中心根据检测统计量T_ξ的表达式，分析主用户信号不存在的情况下的 检测统计量T_ξ的概率分布。

5.1)上三角矩阵L的对角元素的概率分布

在主用户信号不存在的情况下，归一化协方差矩阵R'_x为Wishart矩阵，上三角矩 阵L的对角元素l_ii相互独立，且非负数，其中i＝1,...,M，且服从自由度为N-i+1的 卡方分布，用公式表示为：

$\left(\begin{array}{c}{l}_{\mathrm{ii}}\ge 0\\ l_{\mathrm{ii}}^2~{\chi }_{N-i+1}^2\end{array}\right),$

其中，表示自由度为N-i+1的卡方分布；

5.2)上三角矩阵L的最大特征值的概率分布的累积量分布函数

由上三角矩阵的特性可得，上三角矩阵L的特征值即为L的对角线元素，由步骤 5.1)的分析可得，上三角矩阵L的对角元素为非负数且相互独立，因此，上三角矩阵L 的特征值也非负数且相互独立，

根据以上分析，可得最大特征值ζ₁的概率分布的累积量分布函数用公式表示为：

$F_{{\zeta }_1}\left(y\right)=P(l_{11}\le y,l_{22}\le y,...,l_{\mathrm{MM}}\le y)=\underset{i=1}{\overset{M}{\Pi }}P(l_{\mathrm{ii}}^2\le y^2)=\underset{i=1}{\overset{M}{\Pi }}{\int }_0^{y^2}\frac{x^{\frac{N-i-1}{2}}{e}^{-\frac{x}{2}}}{2^{\frac{N-i+1}{2}}\Gamma \left(\frac{N-i+1}{2}\right)}\mathrm{dx},$

其中，为最大特征值ζ₁的概率分布的累积量分布函数，其中y取值为 (-∞,+∞)，N为每个检测用户的采样点数，M为检测用户数，Γ(·)为伽玛函数。

5.3)上三角矩阵L的最小特征值的概率分布

根据步骤5.2)得到的上三角矩阵L的对角元素即为特征值，以及上三角矩阵L的对 角元素的非负性和独立性，可得最小特征值ζ_M的概率分布的累积量分布函数，用公式 表示为：

$F_{{\zeta }_M}\left(y\right)=1-P(l_{11}^2\ge y^2,l_{22}^2\ge y^2,...,l_{\mathrm{MM}}^2\ge y^2)\text{=1-}\underset{i=1}{\overset{M}{\Pi }}{\int }_{y^2}^{+\infty }\frac{x^{\frac{N-i-1}{2}}{e}^{-\frac{x}{2}}}{2^{\frac{N-i+1}{2}}\Gamma \left(\frac{N-i+1}{2}\right)}\mathrm{dx},$

其中表示最小特征值ζ_M的概率分布的累积量分布函数，其中y取值为 (-∞,+∞)，N为每个检测用户的采样点数，M为检测用户数，Γ(·)为伽玛函数。

5.4)最大特征值ζ₁和最小特征值ζ_M的高斯近似分布

根据中心极限定理，可将最大特征值ζ₁和最小特征值ζ_M的概率分布近似为高斯分 布，用表示最大特征值ζ₁的均值，表示最大特征值ζ₁的方差，表示最小特征 值ζ_M的均值，表示最小特征值ζ_M的方差，其公式分别表示为：

$\left(\begin{array}{c}{u}_{{\zeta }_1}=\underset{-\infty }{\overset{+\infty }{\int }}y{F}_{{\zeta }_1}^{\prime }\left(y\right)\mathrm{dy}\\ {\sigma }_{{\zeta }_1}^2=\underset{-\infty }{\overset{+\infty }{\int }}{(y-u_{{\zeta }_1})}^2{F}_{{\zeta }_1}^{\prime }\left(y\right)\mathrm{dy}\\ u_{{\zeta }_M}=\underset{-\infty }{\overset{+\infty }{\int }}y{F}_{{\zeta }_M}^{\prime }\left(y\right)\mathrm{dy}\\ {\sigma }_{{\zeta }_M}^2=\underset{-\infty }{\overset{+\infty }{\int }}{(y-u_{{\zeta }_M})}^2{F}_{{\zeta }_M}^{\prime }\left(y\right)\mathrm{dy}\end{array}\right),$

其中，为步骤5.2)得到的最大特征值ζ₁的概率分布的累积量分布函数， 为步骤5.3)得到的最小特征值ζ_M的概率分布的累积量分布函数；

5.5)检测统计量T_ξ的概率分布

假设最大特征值ζ₁和最小特征值ζ_M相互独立，根据ζ₁和ζ_M的高斯近似分布，可 以得到检测统计量T_ξ的概率分布为：

$F_{T_{\xi }}\left(y\right)=\Phi \left(\frac{u_{{\zeta }_M}y-u_{{\zeta }_1}}{\sqrt{{\sigma }_{{\zeta }_M}^2{y}^2+{\sigma }_{{\zeta }_1}^2}}\right),$

其中为主用户信号不存在情况下，检测统计量T_ξ的概率分布的累积量分布 函数，Φ(·)为标准正态分布的累积量分布函数，其表示如下：

$\Phi \left(x\right)=\underset{-\infty }{\overset{x}{\int }}\frac{1}{2\pi }{e}^{-\frac{u^2}{2}}\mathrm{du},$

其中，自变量x的取值为(-∞,+∞)。

步骤6，处理中心根据检测统计量T_ξ的概率分布计算判决门限γ_ξ。

根据步骤5得到的检测统计量T_ξ的概率分布，计算检测门限γ_ξ如下：

${\gamma }_{\xi }=\frac{u_{{\zeta }_M}{u}_{{\zeta }_1}}{u_{{\zeta }_M}^2-{\Phi }^{-1}(1-P_{\mathrm{fa}}){\sigma }_{{\zeta }_M}^2}+\frac{{\Phi }^{-1}(1-P_{\mathrm{fa}})\sqrt{{\sigma }_{{\zeta }_1}^2{u}_{{\zeta }_M}^2+{\sigma }_{{\zeta }_M}^2{u}_{{\zeta }_1}^2-{\Phi }^{-1}(1-P_{\mathrm{fa}}){\sigma }_{{\zeta }_1}^2{\sigma }_{{\zeta }_M}^2}}{u_{{\zeta }_M}^2-{\Phi }^{-1}(1-P_{\mathrm{fa}}){\sigma }_{{\zeta }_M}^2},$

其中P_fa表示检测的虚警概率，取值为(0,1)，Φ^-1(·)为标准正态分布的累积量分布函数 Φ(·)的逆函数。

步骤7，处理中心根据检测统计量T_ξ和检测门限γ_ξ，判决主用户信号是否存在。

处理中心通过对步骤4中计算得到的检测统计量T_ξ与步骤6中计算得到的检测门限 γ_ξ进行比较，判决主用户信号是否存在：

当T_ξ≥γ_ξ时，判决为主用户存在，即当前频段频谱已被某用户占用，否则，判决为 主用户不存在，即当前频段频谱为空闲状态，允许检测用户利用。

本发明的频谱检测效果可以通过以下仿真进一步说明：

A、仿真条件

主用户信号为BPSK信号，采用的噪声为均值是0，方差是1的高斯白噪声，仿真 方法为10000000次的蒙特卡洛仿真。对于仿真1-3，检测用户数和采样点数分别设置 为10和20；20和40以及40和100三种情况，虚警概率设置为0.1。对于仿真4，信 噪比设置为从-8dB到8dB，检测用户数和采样点数分别设置为40和100，虚警概率设 置为0.1。对于仿真5，检测用户数和采样点数分别设置为检测用户数和采样点数分别 设置为40和100，性噪比设为0dB。

B、仿真内容

仿真1：在主用户信号不存在情况下，针对理论和仿真，对归一化协方差矩阵的最 大特征值的概率分布的累积量分布进行对比，结果如图2所示，其中，“仿真CDF”表 示本发明的最大特征值的实验累积量分布函数曲线。“近似CDF”表示本发明的最大特 征值的理论累积量分布函数曲线。“10和20”，“20和40”，“40和100”分别表示三种 不同的检测用户数和采样点数的设置组合。

仿真2：在主用户信号不存在情况下，针对理论的和仿真，对归一化协方差矩阵的 最小特征值的概率分布的累积量分布进行对比，结果如图3所示，其中，“仿真CDF” 表示本发明的最大特征值的实验累积量分布函数曲线。“近似CDF”表示本发明的最大 特征值的理论累积量分布函数曲线。“10和20”，“20和40”，“40和100”分别表示三 种不同的检测用户数和采样点数的设置组合。

仿真3：在主用户信号不存在情况下，针对理论的和仿真，对归一化协方差矩阵的 最大特征值和最小特征值之比的概率分布的累积量分布进行对比，结果如图4所示， 其中，“仿真CDF”表示本发明的最大特征值的实验累积量分布函数曲线。“近似CDF” 表示本发明的最大特征值的理论累积量分布函数曲线。“10和20”，“20和40”，“40和 100”分别表示三种不同的检测用户数和采样点数的设置组合。

仿真4：将本发明方法，现有的基于最大最小特征值之比的协作检测方法以及现有 的基于最大特征值和矩阵迹之比的协作检测方法的抗噪声性能进行对比，结果如图5 所示。其中，“最大特征值与迹算法”表示基于最大特征值和矩阵迹之比值的协作检测 方案，“最大最小特征值算法”表示基于最大和最小特征值之比值的协作检测方案，“提 出方法”表示本发明方法。

仿真5：将本发明方法，现有的基于最大最小特征值之比的协作检测方法以及现有 的基于最大特征值和矩阵迹之比的协作检测方法的ROC曲线进行对比，结果如图6所 示。其中，“最大特征值与迹算法”表示基于最大特征值和矩阵迹之比值的协作检测方 案，“最大最小特征值算法”表示基于最大和最小特征值之比值的协作检测方案，“提 出的方案”表示本发明方法。

C、仿真结果

由图2，图3和图4可得，本发明在检测用户数和采样点数较少的情况下，得到的 实验累积量分布曲线和理论累积量分布曲线基本吻合，因此，本发明的检测门限精准 性高，且需要检测用户数和采样点数少，可以在实际中广泛应用。

由图5和图6可见，当性噪比在-8dB到8dB间时，本发明的检测性能好于基于最 大特征值和最小特征值之比值的协作检测方案以及基于最大特征值和矩阵迹之比值的 协作检测方案，这说明本发明在实际中可以得到更广泛的应用。

综合上述仿真结果和分析可得，本发明所需检测用户数和采样点数少，复杂度低， 检测门限精确度高，检测性能比现有的最大特征值和最小特征值之比值的协作检测方 案和基于最大特征值和矩阵迹之比值的协作检测方案的性能好，可在实际中得到更好 地应用。