一种基于HDP‑NSHMM的频谱感知方法与流程

本发明涉及认知无线电频谱感知技术，尤其是涉及一种基于hdp-nshmm的频谱感知方法。

背景技术：

频谱感知是认知无线电中的核心技术，需要实时监测周围的无线环境，为非授权用户给出可用频谱资源，同时保证及时发现授权用户对当前频段的占用，避免造成干扰。因此，频谱感知的准确性对于认知无线电网络起着非常关键的作用。针对频谱感知技术的大量研究工作提供了许多基于信号处理的检测方法，这些检测方法主要可以划分为非协作频谱感知与协作频谱感知两大类。

使用非协作频谱感知的认知无线电设备可以独立选择不同的检测方法，并各自处理本地的感知结果，作出频谱判决。这样的频谱感知方法无需不同用户之间的信息交互，因而更加简单、快捷、便于实现。常用的非协作检测方法主要有能量检测、匹配滤波器检测、循环平稳特征检测等。虽然非协作感知方法无需信息交互，易于实现，然而由于没有考虑空间分集信息，这些检测方法是有局限性的。因此，许多频谱感知算法选择将多个用户的频谱感知结果进行融合的协作检测方法。协作感知能够减少单独的认知无线电设备在频谱感知的检测过程中面对的噪声和信道衰落等带来的不确定性，从而提高频谱检测性能，得到更加准确的判决结果。

根据认知无线电用户数据融合策略的选择，协作感知主要可以分为集中式协作感知以及分布式协作感知两大类。在集中式协作感知中，每一个参与协作感知的认知无线电用户将各自的感知情况经过多跳通信发送到融合中心，由中心基站基于所有用户的感知结果作出统一的频谱判决。然而，由于每个用户都需要与融合中心进行感知数据的传送和判决决策的接收，这种融合方式会带来很大的通信开销，同时，对于数据融合中心的处理能力和处理速度也有较高的要求。为了减少开销以及对融合中心的要求，可以选择每个用户与相邻用户进行局部的信息交换，独立得到频谱判决结果，无需建立融合中心的方法，即分布式协作感知方法。这样的分布式网络成本较低，灵活性也更好。

非平稳隐马尔可夫模型(nonstationaryhmm，nshmm)是一种具有非静态系统隐藏状态转移概率的特殊的隐马尔可夫模型，与传统的hmm不同，nshmm的隐藏状态转移概率是随状态保持时间的变化而改变的。所以分层狄利克雷过程-非平稳隐马尔可夫模型(hdp-nshmm)将更适用于系统隐藏状态转移概率随当前时刻系统状态保持时间变化的应用场景。

技术实现要素：

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于hdp-nshmm的频谱感知方法。

本发明的目的可以通过以下技术方案来实现：

一种基于hdp-nshmm的频谱感知方法，包括以下步骤：

s1、收集所有历史时刻的频谱感知观察数据，做初始化处理；

s2、采用分层狄利克雷过程-非平稳隐马尔可夫模型对历史感知数据进行融合和聚类，在聚类循环中，当信道状态保持时间增长到大于设定时间时，减小自转移偏移参数，使得选择转移到不同的信道状态的可能性增大；

s3、采用伽马分布作为不同信道状态的频谱感知数据指数分布参数的共轭先验来构造贝叶斯模型，利用归为同一信道状态聚类的所有感知数据更新每一类的指数分布参数的伽马分布超参数；

s4、通过更新的伽马分布超参数估计每一类信道状态的功率值，将功率估计值与预设的门限值进行对比，得到频谱判决结果。

优选的，所述步骤s1的初始化处理具体包括：

给每个历史时刻分配初始信道状态类别，统计初始计数和到各历史时刻为止聚类结果保持为初始信道状态的时间序列的长度及保持的起始、结束时刻。

优选的，所述步骤s2中在每一次聚类循环中，对每个时刻t分别进行以下步骤：

s21、获取信道隐藏状态聚类结果及融合参数，对过渡分布β进行抽样，获取时刻t在本次聚类循环前的聚类结果的统计数据；

s22、采样得到t时刻新的信道状态聚类结果，如果选择了新的类别，则更新时刻t的信道状态聚类结果及其相关的统计数据，更新过渡分布，将时刻t对应的频谱感知数据加入相应类的信道状态融合数据中；

s23、将新的信道状态聚类结果与旧的聚类结果比较，若不相同，则更新t前后时刻的时间序列的状态，保持序列长度和起始时间序号；

s24、检查此时是否存在空类别，若存在空类，则移除空类，将总类别数目减一，且将βk加入其中，βk表示选择隐藏状态k的概率，为狄利克雷过程得到的断棍构造方法中的参数，服从βk|α,h～beta(1,α)，表示选择新的隐藏状态的概率；

s25、对融合参数中的辅助变量进行抽样。

优选的，所述步骤s21具体包括：

s211、收集时刻t以外其余时刻的信道隐藏状态聚类结果以及每一类信道状态的融合参数；

s212、对过渡分布β按照狄利克雷分布进行抽样；

s213、将时刻t本次聚类循环前的聚类结果从相关的统计数据中去除，将t时刻的频谱感知数据从该聚类结果的融合数据中移除。

优选的，所述步骤s22具体包括：

s221、分别计算当时刻t选择已有类别或新类别时，到t时刻为止具有连续相同聚类结果的时间序列长度；

s222、对已有的信道状态聚类类别，通过时间序列聚类类别选择的先验分布和似然函数，结合频谱感知数据，计算转移到每个已有信道状态类别的概率和选择新的聚类类别的概率；

s223、采样得到时刻t新的信道状态聚类结果，如果选择了新的类别，则更新时刻t的信道状态聚类结果及其相关的统计数据，更新过渡分布；

s224、将时刻t最新的信道状态聚类结果加入相关的统计数据中，并且将时刻t对应的频谱感知数据加入相应类的信道状态融合数据中。

优选的，所述步骤s222中已有的信道状态聚类类别转移到每个已有信道状态类别的概率为：

其中，yt表示t时刻频谱感知观察数据集合，zt表示每个时刻初始分配信道状态类别，τt表示到t时刻为止聚类结果保持为zt的时间序列的长度，k表示第k个信道状态类别，k表示当前已有的信道类别总数目，κ表示自转移偏移参数，δ(k,zt-1)和δ(k,zt+1)表示冲激函数，和nk′·分别表示t-1、t+1以及所有时刻t＝1,2,...,t从状态类别k′转移出的次数，α和βk为狄利克雷过程得的断棍构造方法中的参数，α为第二层狄利克雷的集中参数，βk表示选择第k个信道状态类别的概率，βk|α,h～beta(1,α)，p(yt；ak,bk)表示似然函数，ak、bk表示伽玛分布的两个超参数；

选择新的聚类类别k+1的概率：

其中，表示t+1时刻聚类结果序号为zt+1的概率，表示选择新的隐藏状态类别的概率。

优选的，所述步骤s223中更新过渡分布具体包括：

对分布进行抽样，令选择新的聚类类别k′的概率为其中表示选择新的隐藏状态类别的概率，βk表示选择第k个信道状态类别的概率，k表示原有的隐藏状态类别数目。

优选的，所述步骤s4具体包括：

通过更新的伽马分布超参数估计每一类信道状态的功率值，将功率估计值与预设的门限值进行对比，如果功率估计值小于门限值，则认为信道可用，如果功率估计值大于门限值，则认为信道被pu占用，不能被su使用；若同时存在多个可用信道，选择接入功率估计值较小的信道。

与现有技术相比，本发明利用hdp-nshmm建立与pu占用信道的状态和持续时间关联的接入模型，通过与聚类类别保持时间相关的自转移偏移参数调整状态自转移概率，对于信道状态的变化能够更准确的判断，通过固定的状态类别可避免冗余状态的出现，提高了历史感知数据聚类准确性，从而具有更高的感知性能，提高了频谱判决准确性。

附图说明

图1为概率图模型hdp-nshmm的示意图；

图2为实施例二中马尔可夫模型下主用户对信道占用情况；

图3为实施例二中非平稳马尔可夫模型下主用户对信道占用情况；

图4为实施例二中不固定类别k时初始化的信道隐藏状态分配结果；

图5为实施例二中不固定类别k时1次循环后的信道隐藏状态分配结果；

图6为实施例二中不固定类别k时10次循环后的信道隐藏状态分配结果；

图7为实施例二中不固定类别k时100次循环后的信道隐藏状态分配结果；

图8为实施例二中固定类别k时初始化的信道隐藏状态分配结果；

图9为实施例二中固定类别k时1次循环后的信道隐藏状态分配结果；

图10为实施例二中固定类别k时10次循环后的信道隐藏状态分配结果；

图11为实施例二中固定类别k时100次循环后的信道隐藏状态分配结果；

图12为实施例二中每类信道状态对应功率估计值；

图13为实施例二中实际的信道可用情况；

图14为实施例二中信道状态判决结果；

图15为实施例二中不同频谱感知方法在接收机信噪比为0db时的roc对比曲线；

图16为实施例二中不同频谱感知方法在接收机信噪比为05db时的roc对比曲线；

图17为实施例二中不同频谱感知方法在接收机信噪比为10db时的roc对比曲线；

图18为实施例二中不同频谱感知方法在虚警概率为0.2时的正确检测概率。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

实施例一

本方法将分层狄利克雷过程与非平稳隐马尔可夫模型相结合，实现历史感知数据的自动聚类。

图1所示为本方法hdp-nshmm的概率图模型，将常量κ改变为随时间变化的函数变量κ(τ)使得状态转移概率与状态持续时间相关，成为非时间平稳的模型。矩形方框中的k表示框中的所有τt-1个隐藏状态类别序号都为k，即满足本申请提出一种基于hdp-nshmm的频谱感知方法，包括以下步骤：

步骤1、数据初始化过程

步骤(11)、收集所有t个历史时刻的频谱感知观察数据集合{yt|t＝1,2,...,t}；

步骤(12)、在大规模认知无线电网络频谱感知的场景中，首先根据用户的空间分布将当前可能具有相同信道状态分布参数的用户分为同一组，具体方法为将具有相近指数分布参数λk的历史感知信号分为一组，并在组内进行分布式的频谱感知数据融合，为每个时刻初始分配信道状态类别{zt|t＝1,2,...,t}；

步骤(13)、统计初始计数nij，即所有时刻t＝1,2,...,t从状态类别i转移到类别j的次数，统计计数ni·，即所有时刻t＝1,2,...,t从状态类别i转移出的次数。其中，i,j＝1,2,...,k，k表示初始化的隐藏状态类别数目；

步骤(14)、统计τ1,τ2,…,τt，其中，τt表示到t时刻为止，聚类结果保持为zt的时间序列长度；

步骤(15)、记录每个时刻所在持续状态的起始时刻ts1,ts2,…,tst以及结束时刻te1,te2,…,tet，其中，tst和tet分别表示t时刻所在的具有连续相同聚类结果的时间序列的起始时间和结束时间。

步骤2、历史感知数据的融合与聚类

在当前状态的保持时间较短时，设置较大的自转移偏移参数，以保证其状态不会随时间快速变化，产生冗余类别，而随着当前状态保持时间的增长，减小自转移偏移参数，从而减小状态的自转移概率，使其更有可能选择转移到不同的状态，最终提高频谱判决结果的准确性，在每一次聚类循环中，对每个时刻t分别进行以下步骤：

步骤(21)、获取信道隐藏状态聚类结果及融合参数，对过渡分布进行抽样，获取时刻t在本次聚类循环前的聚类结果的统计数据；

步骤(211)、收集时刻t以外其余时刻的信道隐藏状态聚类结果{zt|t＝1,2,...,t-1,t+1...,t}以及每一类信道状态的融合参数{m·k,ak,bk|k＝1,2,...,k}；k表示第k个信道状态类别。

步骤(212)、对过渡分布β按照狄利克雷分布进行抽样：

β～dir(m·1,m·2,...,m·k,γ)

步骤(213)、假设t时刻本次聚类循环开始之前的聚类结果为z′t＝k′，将其从k′类相关的统计数据中去除，即将计数和nk′·减一，和nk′·分别表示t-1、t+1以及所有时刻t＝1,2,...,t从状态类别k′转移出的次数，并将t时刻的频谱感知数据从k′类的融合数据中移除，同时记录下旧的聚类结果k′；

步骤(22)、采样得到t时刻新的信道状态聚类结果，如果选择了新的类别，则更新时刻t的信道状态聚类结果及其相关的统计数据，更新过渡分布，将时刻t对应的频谱感知数据加入相应类的信道状态融合数据中，具体包括：

步骤(221)、分别计算当t时刻选择k个已有类别和新类别时，到t时刻为止具有连续相同聚类结果zt的时间序列长度τt(k)；

步骤(222)、对于k个当前已有的信道状态聚类类别，时间序列聚类类别选择的先验分布公式为：

似然函数的公式为：

其中，λ＝{λk|wlt/π≤k≤wht/π}为信道的状态参数集合，p(y/λ)表示y在λ已知情况下的条件概率，lomax(y(k)；ak,bk)表示y(k)在超参数ak和bk条件下lomax分布，cg(y(k)；l＝1,ak,bk)表示复合伽玛分布在l＝1时的一种特别的情况，该分部即为观察数据y的似然函数，wl表示低频信号带宽，wh表示高频信号带宽，l表示gamma分布的特定参数。

结合频谱感知数据，计算转移到每一个已有信道状态类别k的概率：

其中，yt表示t时刻频谱感知观察数据集合，zt表示每个时刻初始分配信道状态类别，τt表示到t时刻为止聚类结果保持为zt的时间序列的长度，k表示第k个信道状态类别，k表示当前已有的信道类别总数目，κ表示自转移偏移参数，δ(k,zt-1)和δ(k,zt+1)表示冲激函数，和nk′·分别表示t-1、t+1以及所有时刻t＝1,2,...,t从状态类别k′转移出的次数，α和βk为狄利克雷过程得的断棍构造方法中的两个参数，α为第二层狄利克雷的集中参数，βk|α,h～beta(1,α)，p(yt；ak,bk)表示似然函数，ak、bk表示伽玛分布的两个超参数；

计算选择新的聚类类别k+1的概率：

其中，表示t+1时刻聚类结果序号为zt+1的概率，表示除现有隐藏状态外新的隐藏状态类别的概率；

特别地，对于初始时间的感知数据，即t＝1时，有：

对于时间序列的最后一个感知数据，即t＝t时，有：

步骤(223)、采样得到时刻t新的信道状态聚类结果zt：

如果选择了新的类别，增加k值，更新过渡分布β：对分布进行抽样，令选择新的聚类类别k′的概率为其中表示选择新的隐藏状态类别的概率，βk表示选择第k个信道状态类别的概率；

步骤(224)、将时刻t最新的信道状态聚类结果加入相关的统计数据中，即将计数和nk′·加一，并且将时刻t对应的频谱感知数据加入相应类的信道状态融合数据中；

步骤(23)、将新的信道状态聚类结果与旧的聚类结果比较，若不相同，则从tst-1到tet+1的时间序列的状态，保持序列长度和起始时间序号；

步骤(24)、检查此时是否存在空类别，若存在空类即对所有k＝1,2,...,k，检查是否存在nk·＝0且n·k＝0，若存在空类，则移除空类k，将总类别数目k减一，且将βk加入以提高方法执行效率；

步骤(25)、对辅助变量m进行抽样。即对于步骤(12)分组q中的每一个信道状态k＝1,2,...,k，初始化mqk＝0，表示q组中信道状态k的个数为0，为分组q每一个信道状态为k的观察值yt，从下列分布中抽样：

其中，n表示信道状态的总个数，从0到njk-1逐个增加n，设x为从分布中取出的抽样值，如果x＝1，则mqk加1。

步骤3、历史感知数据分布参数更新

采用伽马分布作为不同信道状态的频谱感知数据指数分布参数λk的共轭先验来构造贝叶斯模型：

γ(ak)表示以ak为参数的gamma函数；

上式中，表示n个抽样后的频域信号y(k)的观察值的平均，p(λk|{y(k)}n)表示在取得了y(k)的n个观察值之后λk的后验分布，exponential({y(k)}n|λk)表示指数分布，gamma(λk；ak,bk)表示伽马分布，a′k、b′k表示ak和bk更新后的值。

利用归为同一信道状态聚类的所有感知数据更新每一类的指数分布参数的伽马分布超参数(ak,bk)。

上式中，为可以衡量信道情况的的估计值，e[(1/λk)|{y(k)}n]表示在得到n个观察值y(k)的情况下，1/λk的期望值。

步骤4、频谱状态判决

步骤(41)、信道状态聚类以及聚为同一类的频谱感知数据融合结束之后，可以得到每一类信道状态参数λk更新后的分布参数(ak,bk)，可对每一类信道状态的功率值进行估计：

步骤(42)、将功率估计值与预设的门限值进行对比，如果功率估计值小于门限值，则认为信道可用，如果功率估计值大于门限值，则认为信道被pu占用，不能被次用户(su)使用；若同时存在多个可用信道，选择接入功率估计值较小的信道。

实施例二

本实施例验证在主用户(pu)对信道占用与释放情况的转移概率随状态持续时间变化的场景下，本申请提出的频谱感知方法的聚类效果和频谱判决性能。

仍考虑小型认知无线电网络的场景，pu的信号经过瑞利信道，并按自由空间传播模型衰减，不同的是，pu对信道的占用和释放情况服从非平稳的马尔可夫模型，其初始的状态转移概率仍设置为p0(0/0)＝0.975，p0(1/0)＝0.025；p0(0/1)＝0.05，p0(1/1)＝0.95，而当前时刻状态已保持时间大于预设值时，将减少状态自转移概率。在本实施例的仿真设置中，令预设保持时间为10个时间序列长度，保持时间大于预设时间后，按与初始转移概率相同的比例减小状态自转移概率，即令pt(0/0)＝pt-1(0/0)·p0(0/0)或pt(1/1)＝pt-1(1/1)·p0(1/1)。

基于这一状态转移概率随持续时间变化的非平稳马尔可夫模型与固定转移概率的马尔可夫模型产生的pu对信道的占用与释放情况分别如图2和图3所示。横坐标为时间序列序号，总长度设置为t＝200，纵坐标为pu对信道的占用情况，“1”表示占用信道，“0”表示释放信道。从图中可以看出，与固定的状态转移概率相比，非平稳的状态转移概率会使得pu对信道的占用和释放持续时间不至过长，且更具有规律性。

下面验证非平稳马尔可夫模型场景产生的频谱感知数据的聚类效果。对时间序列长度t＝200的频谱感知数据进行聚类，su接收机信噪比设置为10db。聚类算法的参数初始化为：α＝2，γ＝2，α和γ是dirichlet分布的两个特定参数，初始状态数目k＝4，a＝l＝2，a、l、b是gamma分布的三个特定参数，表示所有时刻的频谱感知观察数据的平均值状态类别数目k可变。考虑单用户对两个信道的感知场景，两个pu对两个信道的占用情况都符合上述相同非平稳马尔可夫模型。频谱感知数据的隐藏状态聚类情况如图4～7所示，依次为初始化、1次循环后、10次循环后以及100次循环后的信道隐藏状态分配结果。由于状态自转移偏移参数κ(τ)将随着状态持续时间而减小，使得选择新类别的概率增加，从而更容易产生冗余状态，影响聚类的准确性。因此，在已知隐藏信道状态类别的前提下，可以固定聚类类别数目k，计算每个时刻的感知数据属于每一个已有类别k的概率，并从k个已有类别中进行抽样，确定聚类结果。

对于两个pu的仿真场景，将类别数目固定为k＝4，此时初始化、1次循环后、10次循环后以及100次循环后的信道隐藏状态分配结果如图8～11所示。循环结束后，计算出的每一类的功率估计值如图12所示，横坐标表示聚类类别，纵坐标表示功率估计值sk。在图12中，类别1表示第一个信道可用，第二个信道不可用，记为(0,1)；类别2表示两个信道都可用，记为(0,0)；类别3表示两个信道都不可用，记为(1,1)；类别4表示第一个信道不可用，第二个信道可用，记为(1,0)。将两个信道的四种可用情况(0,0)、(0,1)、(1,0)、(1,1)分别记为信道状态1、状态2、状态3、状态4。信道状态判决结果与实际的信道可用情况分别如图13和14所示。

将信道功率估计值与不同的预设门限值比较，可画出正确检测概率与虚警概率形成的roc曲线。对比本申请提出的基于hdp-nshmm的固定类别聚类与频谱判决方法与固定自转移偏移参数的stickyhdp-hmm方法(专利“一种基于多用户历史感知数据挖掘的协作频谱感知方法”，申请公布号cn106972899a)、不加入自转移偏移参数的hdp-hmm方法以及能量检测、匹配滤波器检测、循环平稳特征检测在相同接收机信噪比下的方法性能。图15～17分别为接收机信噪比为0db、5db、10db时的roc曲线对比。图18为虚警概率为0.2处不同信噪比下各个频谱感知方法的正确检测概率对比图。

在信噪比为0db、5db、10db时，本方法在虚警概率为0.2处的正确检测概率比能量检测分别提高了52％、31％、9％。本方法在较高信噪比下，虚警概率为0.2时的正确检测概率低于匹配滤波器检测和循环平稳特征检测，而匹配滤波器检测需要预先获知pu传输信号的先验知识，循环平稳特征检测系统复杂度较高，且需要较长的检测时间，在实际应用中具有一定的局限性。在较低信噪比下，由于利用历史数据对具有较大不确定性的频谱感知数据进行了聚类融合，本方法在虚警概率为0.2处具有最高的正确检测概率。在信噪比为0db时，本方法在虚警概率为0.2处的正确检测概率比匹配滤波器检测和循环平稳特征检测分别提高了15％和31％。而stickyhdp-hmm方法在聚类过程中使用了固定的自转移偏移参数，hdp-hmm算法不加入自转移偏移参数，无法区分状态自转移概率，都不适用于状态转移概率随持续时间变化的应用场景，在低信噪比下容易导致错误的聚类结果，从而进行错误的数据融合，在相同的虚警概率下具有较低的正确检测概率。本方法由于利用了灵活的自转移偏移参数，对于信道状态的变化有着更准确的判断，且固定的状态类别可避免冗余状态的出现，提高了历史感知数据聚类准确性，从而具有更高的感知性能。在信噪比为0db、5db、10db时，本方法在虚警概率为0.2处的正确检测概率比stickyhdp-hmm方法分别提高了63％、18％、4％。