CN112202513B - 基于生灭过程和粘性隐马尔可夫模型的多带频谱感知方法 - Google Patents

Abstract

本发明公开了一种基于生灭过程和粘性隐马尔可夫模型的多带频谱感知方法，其计算每个频带对应的接收信号功率，并将其作为隐马尔可夫模型的观测数据，再确定观测数据的隐藏状态、属于各类隐藏状态的概率，计算状态转移概率矩阵；隐马尔可夫模型中引入粘性因子得到粘性隐马尔可夫模型；在每次迭代下计算所有观测数据的隐藏状态的聚类结果，通过生灭过程更新聚类结果，并计算隐藏状态的类别数，计算状态转移概率矩阵及均值和精度；迭代结束后根据各类隐藏状态的功率估计值和观测数据的隐藏状态在最后一次迭代下的值计算频带的功率估计值，通过与门限值比较，确定频带是否被授权用户占用；优点是对多个频带同时进行频谱感知不受频带稀疏性的影响。

Description

基于生灭过程和粘性隐马尔可夫模型的多带频谱感知方法

技术领域

本发明涉及一种认知无线电中的频谱感知技术，尤其是涉及一种基于生灭过程和粘性隐马尔可夫模型的多带频谱感知方法。

背景技术

与第四代移动通信技术相比，第五代移动通信(5G)技术可将数据速率提高到10Gbit/s、将延迟降低到1毫秒、连接设备数量增加100倍。实现这些需求依赖大量的频谱资源，但是可用的频谱资源是有限的，且基本上已经被分配完了。为了解决无线网络频谱资源短缺的问题，目前业界的常见思路是引入认知无线电技术来提高频谱资源的利用率。与传统的频带被单一用户授权占用的系统不同，无线网络通过认知无线电频谱感知技术可以从环境中智能地检测频带占用情况，从而使认知用户能智能接入空闲的授权频带。

在实际应用中，为了充分提高所有频带的频谱利用率，多频带频谱感知方案越来越多的得到了人们的关注。多频带频谱感知方案对多个频带进行频谱感知，当授权用户信号再次出现时，认知用户能快速切换到其它空闲频带。但是，现有的多频带频谱感知方案经常利用频带稀疏性来实现检测，当频带不具有稀疏性时具有较差的性能。

发明内容

本发明所要解决的技术问题是提供一种基于生灭过程和粘性隐马尔可夫模型的多带频谱感知方法，其不管频带是否具有稀疏性，均能够对多个频带同时进行频谱感知，且感知性能好。

本发明解决上述技术问题所采用的技术方案为：一种基于生灭过程和粘性隐马尔可夫模型的多带频谱感知方法，其特征在于具体包括以下步骤：

步骤一：在认知无线电系统中，对连续的L个频带内的信号进行采样，且对每个频带内的信号进行等时间间隔地采样，共采样得到N个样本，将对第j个频带内的信号进行采样得到的第n个样本记为r_j(n)；然后计算每个频带对应的接收信号功率，将第j个频带对应的接收信号功率记为x_j，即为对第j个频带内的信号进行采样得到的所有样本的平均功率，

其中，L、N、j和n均为正整数，L＞1，100≤N≤1000，1≤j≤L，1≤n≤N，符号“||”为求模符号，x_j服从高斯分布，即

表示噪声功率，

表示第j个频带内授权用户的信号功率，

表示第j个频带未被授权用户占用，

表示第j个频带已被授权用户占用，

表示x_j服从均值为

方差为

的高斯分布，

表示x_j服从均值为

方差为

的高斯分布；

步骤二：将每个频带对应的接收信号功率作为隐马尔可夫模型中的观测数据，即隐马尔可夫模型中的第j个观测数据为x_j；然后确定隐马尔可夫模型中的每个观测数据所对应的一个隐藏状态，将x_j对应的隐藏状态记为z_j，z_j的取值为区间[1,K]中的一个值，若z_j的取值为1则认为x_j属于第1类隐藏状态，若z_j的取值为k则认为x_j属于第k类隐藏状态，若z_j的取值为K则认为x_j属于第K类隐藏状态；接着计算隐马尔可夫模型中的每个观测数据属于各类隐藏状态的概率，将x_j属于第k类隐藏状态的概率记为

最后计算隐马尔可夫模型中的状态转移概率矩阵，记为Q，

其中，K和k均为正整数，K表示隐马尔可夫模型中的隐藏状态的类别数，2≤K≤10，1≤k≤K，

表示x_j服从的高斯分布的概率密度函数，其变量为x_j、均值为μ_k、方差为

μ_k表示属于第k类隐藏状态的高斯分布的均值，τ_k表示属于第k类隐藏状态的高斯分布的精度即方差的倒数，Q_1,1、Q_1,2、Q_1,k'、Q_1,K对应表示Q中的第1行第1列的元素、第1行第2列的元素、第1行第k'列的元素、第1行第K列的元素，Q_2,1、Q_2,2、Q_2,k'、Q_2,K对应表示Q中的第2行第1列的元素、第2行第2列的元素、第2行第k'列的元素、第2行第K列的元素，Q_k,1、Q_k,2、Q_k,k'、Q_k,K对应表示Q中的第k行第1列的元素、第k行第2列的元素、第k行第k'列的元素、第k行第K列的元素，Q_K,1、Q_K,2、Q_K,k'、Q_K,K对应表示Q中的第K行第1列的元素、第K行第2列的元素、第K行第k'列的元素、第K行第K列的元素，1≤k'≤K，Q_k,k'表示z_j'-1＝k的条件下z_j'＝k'的概率，2≤j'≤L，z_j'-1表示隐马尔可夫模型中的第j'-1个观测数据x_j'-1对应的隐藏状态，z_j'表示隐马尔可夫模型中的第j'个观测数据x_j'对应的隐藏状态；

步骤三：在隐马尔可夫模型中引入粘性因子，得到粘性隐马尔可夫模型；在粘性隐马尔可夫模型中，初始化属于每类隐藏状态的高斯分布的均值和精度，将μ_k的初始化值记为

将τ_k的初始化值记为

初始化隐藏状态的类别数K，将K的初始化值记为K⁽⁰⁾，K⁽⁰⁾为区间[2,10]内的任意一个正整数；初始化状态转移概率矩阵Q，将Q的初始化值记为Q⁽⁰⁾，Q⁽⁰⁾中的每行中的所有元素的共轭先验分布服从狄利克雷分布，Q⁽⁰⁾中的第k⁽⁰⁾行中的所有元素的共轭先验分布服从的狄利克雷分布为：

其中，

表示Q⁽⁰⁾中的第k⁽⁰⁾行中的所有元素，Dir()表示狄利克雷分布，γ表示狄利克雷分布的参数，κ表示粘性因子，δ(k⁽⁰⁾,1)表示两个参数分别为k⁽⁰⁾和1的克罗内克函数，δ(k⁽⁰⁾,k'⁽⁰⁾)表示两个参数分别为k⁽⁰⁾和k'⁽⁰⁾的克罗内克函数，δ(k⁽⁰⁾,K⁽⁰⁾)表示两个参数分别为k⁽⁰⁾和K⁽⁰⁾的克罗内克函数，

γ+κδ(k⁽⁰⁾,1)表示共轭先验分布服从的狄利克雷分布的第1个元素，γ+κδ(k⁽⁰⁾,k'⁽⁰⁾)表示共轭先验分布服从的狄利克雷分布的第k'⁽⁰⁾个元素，γ+κδ(k⁽⁰⁾,K⁽⁰⁾)表示共轭先验分布服从的狄利克雷分布的第K⁽⁰⁾个元素，1≤k(⁰)≤K(⁰)，1≤k'(⁰)≤K(⁰)；

步骤四：令t表示迭代次数，t的初始值为1；令T表示设定的最大迭代次数，T≥3；

步骤五：计算在第t次迭代下粘性隐马尔可夫模型中的所有观测数据对应的隐藏状态的聚类结果，记为z^(t)，

其中，

表示求使得p(z|x,Q^(t-1),μ^(t-1),τ^(t-1))取最大值时变量z的值，z为粘性隐马尔可夫模型中的所有观测数据对应的隐藏状态构成的向量，z＝[z₁,z₂,…,z_j,…,z_L]，z₁表示第1个观测数据x₁对应的隐藏状态，z₂表示第2个观测数据x₂对应的隐藏状态，z_L表示第L个观测数据x_L对应的隐藏状态，x表示粘性隐马尔可夫模型中的所有观测数据构成的向量，x＝[x₁,x₂,…,x_j,…,x_L]，t＝1时Q^(t-1)即为Q⁽⁰⁾，t≠1时Q^(t-1)表示在第t-1次迭代下粘性隐马尔可夫模型中的状态转移概率矩阵Q的值，t＝1时μ^(t-1)即为μ的初始值μ⁽⁰⁾，μ＝[μ₁,μ₂,…,μ_K]，μ₁表示属于第1类隐藏状态的高斯分布的均值，μ_K表示属于第K类隐藏状态的高斯分布的均值，

表示μ₁的初始化值，

表示μ₂的初始化值，

表示属于第K⁽⁰⁾类隐藏状态的高斯分布的均值

的初始化值，t≠1时μ^(t-1)表示μ在第t-1次迭代下的值，

表示μ₁在第t-1次迭代下的值，

表示μ₂在第t-1次迭代下的值，

表示属于第K^(t-1)类隐藏状态的高斯分布的均值

在第t-1次迭代下的值，t＝1时τ^(t-1)即为τ的初始值τ⁽⁰⁾，τ＝[τ₁,τ₂,…,τ_K]，τ₁表示属于第1类隐藏状态的高斯分布的精度，τ_K表示属于第K类隐藏状态的高斯分布的精度，

表示τ₁的初始化值，

表示τ₂的初始化值，

表示属于第K⁽⁰⁾类隐藏状态的高斯分布的精度

的初始化值，t≠1时τ^(t-1)表示τ在第t-1次迭代下的值，

表示τ₁在第t-1次迭代下的值，

表示τ₂在第t-1次迭代下的值，

表示属于第K^(t-1)类隐藏状态的高斯分布的精度

在第t-1次迭代下的值，t＝1时K^(t-1)即为K的初始化值K⁽⁰⁾，t≠1时K^(t-1)表示K在第t-1次迭代下的值，

表示z的后验概率，根据贝叶斯定理得到

p(z_j|x,Q^(t-1),μ^(t-1),τ^(t-1))表示z_j的后验概率，符号“∝”表示正比，x_j+1表示第j+1个观测数据，x_j+2表示第j+2个观测数据，p(z_j,x₁,x₂,...,x_j|Q^(t-1),μ^(t-1),τ^(t-1))表示z_j,x₁,x₂,...,x_j的联合概率，p(x_j+1,x_j+2,...,x_L|z_j,Q^(t-1),μ^(t-1),τ^(t-1))表示z_j的条件下x_j+1,x_j+2,...,x_L的联合概率，p(z_j,x₁,x₂,...,x_j|Q^(t-1),μ^(t-1),τ^(t-1))和p(x_j+1,x_j+2,...,x_L|z_j,Q^(t-1),μ^(t-1),τ^(t-1))通过前向后向算法计算得到；

步骤六：通过生灭过程更新z^(t)，进而计算在第t次迭代下粘性隐马尔可夫模型中的隐藏状态的类别数K的值K^(t)，具体过程为：

1)统计z^(t)中值等于区间[1,K^(t-1)]中的每个值的元素的总个数，将z^(t)中值等于k^(t-1)的元素的总个数记为

2)按1,…,k^(t-1),…,K^(t-1)从小到大的顺序排列Num₁至

得到统计个数排列序列；

3)若统计个数排列序列中只有一个0值且K^(t-1)≠2，则假设该0值对应区间[1,K^{(t -1)}]中的k^(t-1)，那么将z^(t)中值分别等于(k+1)^(t-1)至K^(t-1)的所有元素的值均减1，将更新后的z^(t)重新记为z^*(t)，并令K^(t)＝K^(t-1)-1，再执行步骤七；若统计个数排列序列中只有一个0值且K^(t-1)＝2，或者统计个数排列序列中有多个0值，则执行步骤4)；若统计个数排列序列中没有0值且K^(t-1)＜10，则从1到L中随机产生一个起始值，记为L_min，从L_min到L中随机产生一个终止值，记为L_max，将z^(t)中的

的值均设为K^(t-1)+1，将更新后的z^(t)重新记为z^*(t)，并令K^(t)＝K^(t-1)+1，再执行步骤七；如果统计个数排列序列为除上述四种情况外的其余情况，则保持z^(t)不变，将z^(t)重新记为z^*(t)，并令K^(t)＝K^(t-1)，再执行步骤七；

4)当统计个数排列序列中的第ω个值为L而第1个至第ω-1个以及第ω+1个至第K^(t-1)个值均为0即z^(t)中的所有元素的值等于ω时，从1到L中随机产生一个起始值，记为L'_min，从L'_min到L中随机产生一个终止值，记为L'_max，将z^(t)中的

的值均设为2，将z^(t)中的

和

的值均设为1，将更新后的z^(t)重新记为z^*(t)，并令K^(t)＝2，再执行步骤七；当统计个数排列序列中有多个0值且任一个非0值不为L时，执行步骤5)；

5)设0值有ξ个，针对第1个0值和第2个0值，假设第1个0值对应区间[1,K^(t-1)]中的k^(t-1)、第2个0值对应区间[1,K^(t-1)]中的(k+υ)^(t-1)，那么将z^(t)中值分别等于(k+1)^(t-1)至(k+υ-1)^(t-1)的所有元素的值均减1；针对第2个0值和第3个0值，假设第2个0值对应区间[1,K^{(t -1)}]中的(k+υ)^(t-1)、第3个0值对应区间[1,K^(t-1)]中的

那么将z^(t)中值分别等于(k+υ+1)^(t-1)至

的所有元素的值均减2；依次类推，针对第ξ-1个0值和第ξ个0值，假设第ξ-1个0值对应区间[1,K^(t-1)]中的(k+θ)^(t-1)、第ξ个0值对应区间[1,K^(t-1)]中的(k+ρ)^(t-1)，那么将z^(t)中值分别等于(k+θ+1)^(t-1)至(k+ρ-1)^(t-1)的所有元素的值均减ξ-1；而后将z^(t)中值分别等于(k+ρ+1)^(t-1)至K^(t-1)的所有元素的值均减ξ；将更新后的z^(t)重新记为z^*(t)，并令K^(t)＝K^(t-1)-ξ，再执行步骤七；

上述，k^(t-1)为区间[1,K^(t-1)]中的第k^(t-1)个值，Num₁表示z^(t)中值等于1的元素的总个数，

表示z^(t)中值等于K^(t-1)的元素的总个数，1≤ω≤K^(t-1)，1≤L_min≤L_max≤L，

对应表示z^(t)中的第L_min个元素、第L_min+1个元素、…、第L_max个元素，1≤L'_min≤L'_max≤L，

对应表示z^(t)中的第L'_min个元素、第L'_min+1个元素、…、第L'_max个元素，

对应表示z^(t)中的第1个元素、第2个元素、…、第L'_min-1个元素，

对应表示z^(t)中的第L'_max+1个元素、第L'_max+2个元素、…、第L个元素，1＜ξ＜K^(t-1)，

K^(t)＝K^(t-1)中的“＝”为赋值符号，

表示z₁在第t次迭代下经过生灭过程后的值，

表示z₂在第t次迭代下经过生灭过程后的值，

表示z_j在第t次迭代下经过生灭过程后的值，

表示z_L在第t次迭代下经过生灭过程后的值；

步骤七：计算在第t次迭代下粘性隐马尔可夫模型中的状态转移概率矩阵Q的值，记为Q^(t)，Q^(t)中的第k^(t)行中的所有元素的共轭先验分布服从的狄利克雷分布为：

Q^(t)中的第k^(t)行中的所有元素的后验分布服从的狄利克雷分布为：

其中，1≤k^(t)≤K^(t)，1≤k'^(t)≤K^(t)，

表示Q^(t)中的第k^(t)行中的所有元素，

表示在第t次迭代下经过生灭过程后从第k^(t)类隐藏状态转移到第1类隐藏状态的观测数据的数量，

表示在第t次迭代下经过生灭过程后从第k^(t)类隐藏状态转移到第k'^(t)类隐藏状态的观测数据的数量，

表示在第t次迭代下经过生灭过程后从第k^(t)类隐藏状态转移到第K^(t)类隐藏状态的观测数据的数量，δ(k^(t),1)表示两个参数分别为k^(t)和1的克罗内克函数，δ(k^(t),k'^(t))表示两个参数分别为k^(t)和k'^(t)的克罗内克函数，δ(k^(t),K^(t))表示两个参数分别为k^(t)和K^(t)的克罗内克函数，

表示后验分布服从的狄利克雷分布的第1个元素，

表示后验分布服从的狄利克雷分布的第k'^(t)个元素，

表示后验分布服从的狄利克雷分布的第K^(t)个元素；

步骤八：利用属于同一类隐藏状态的所有观测数据，根据贝叶斯定理，计算在x和z^*(t)确定后μ在第t次迭代下的值μ^(t)和τ在第t次迭代下的值τ^(t)的后验概率，记为p(μ^(t),τ^(t)|x,z^*(t))，

其中，μ^(t)表示μ在第t次迭代下的值，

表示μ₁在第t次迭代下的值，

表示属于第k^(t)类隐藏状态的高斯分布的均值

在第t次迭代下的值，

表示属于第K^(t)类隐藏状态的高斯分布的均值

在第t次迭代下的值，τ^(t)表示τ在第t次迭代下的值，

表示τ₁在第t次迭代下的值，

表示属于第k^(t)类隐藏状态的高斯分布的精度

在第t次迭代下的值，

表示属于第K^(t)类隐藏状态的高斯分布的精度

在第t次迭代下的值，

表示在第t次迭代下经过生灭过程后属于第k^(t)类隐藏状态的观测数据的数量，

表示在第t次迭代下经过生灭过程后属于第k^(t)类隐藏状态的所有观测数据的平均值，

表示在第t次迭代下经过生灭过程后属于第k^(t)类隐藏状态的第

个观测数据，

表示

服从的高斯分布的概率密度函数，其变量为

均值为

方差为

表示

服从的伽马分布的概率密度函数，其变量为

形状参数为

速率参数为

η₀、m₀、a₀、b₀均为常数；

步骤九：判断t＜T是否成立，如果成立，则令t＝t+1，然后返回步骤五继续迭代；如果不成立，则执行步骤十；其中，t＝t+1中的“＝”为赋值符号；

步骤十：将μ^(t)中的每个元素的值作为对应一类隐藏状态的功率估计值，即将μ^(t)中的第k^(t)个元素的值作为第k^(t)类隐藏状态的功率估计值；然后根据每类隐藏状态的功率估计值和每个观测数据对应的隐藏状态在第t次迭代下经过生灭过程后的值，计算每个频带的功率估计值，将第j个频带的功率估计值记为β_j，若

则β_j等于第k^(t)类隐藏状态的功率估计值；再对每个频带的功率估计值与门限值进行比较，对于β_j，如果β_j小于门限值，则认为第j个频带未被授权用户占用，并将第j个频带作为可用频带；如果β_j大于或等于门限值，则认为第j个频带已被授权用户占用，不能被认知用户使用；其中，门限值是根据给定的虚警概率计算得到的。

与现有技术相比，本发明的优点在于：

1)本发明方法适用于对多个频带内的信号进行多频带频谱感知，在对多频带内的信号进行频谱感知的过程中，本发明方法通过在相邻频带之间增加粘性因子建立粘性隐马尔可夫模型，来增加相邻频带之间的相关性，由于利用到多频带中相邻频带之间的相关性进行多频带频谱检测，因此在相同条件下本发明方法能够充分提高所有频带的频谱利用率。

2)本发明方法由于利用到多个频带中相邻频带之间的相关性，从而不管频带是否具有稀疏性，均能够对多个频带同时进行频谱感知，并具有良好的检测性能。

3)本发明方法由于利用到粘性隐马尔可夫模型特有的对观测数据的有效聚类能力，因此可以实现快速收敛，从而具有较低的计算复杂度。

4)本发明方法由于利用到生灭过程，因此可以实现自动调整频谱状态的类别总数，从而使算法更加灵活。

附图说明

图1为本发明方法的总体实现框图；

图2为信噪比为-10dB，频带占用率为20％(频带具有稀疏性)时本发明方法和使用无粘性的隐马尔可夫模型的方法的平均功率估计误差随迭代次数的变化曲线；

图3为信噪比为-10dB，频带占用率为20％(频带具有稀疏性)时本发明方法和使用无粘性的隐马尔可夫模型的方法的ROC曲线；

图4为信噪比为-10dB，频带占用率为80％(频带不具有稀疏性)时本发明方法和使用无粘性的隐马尔可夫模型的方法的ROC曲线。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

本发明提出的一种基于生灭过程和粘性隐马尔可夫模型的多带频谱感知方法，其总体实现框图如图1所示，其具体包括以下步骤：

步骤一：在认知无线电系统中，对连续的L个频带内的信号进行采样，且对每个频带内的信号进行等时间间隔地采样，共采样得到N个样本，将对第j个频带内的信号进行采样得到的第n个样本记为r_j(n)；然后计算每个频带对应的接收信号功率，将第j个频带对应的接收信号功率记为x_j，即为对第j个频带内的信号进行采样得到的所有样本的平均功率，

其中，L、N、j和n均为正整数，L＞1，100≤N≤1000，1≤j≤L，1≤n≤N，符号“||”为求模符号，N的取值太大(N＞1000)时采样的样本数量太多，会导致运算速度下降，因此只需确保N的取值充分大即可，当N充分大时，根据中心极限定理，x_j服从高斯分布，即

表示噪声功率，

表示第j个频带内授权用户的信号功率，

表示第j个频带未被授权用户占用，

表示第j个频带已被授权用户占用，

表示x_j服从均值为

方差为

的高斯分布，

表示x_j服从均值为

方差为

的高斯分布。

步骤二：将每个频带对应的接收信号功率作为隐马尔可夫模型中的观测数据，即隐马尔可夫模型中的第j个观测数据为x_j；然后确定隐马尔可夫模型中的每个观测数据所对应的一个隐藏状态，将x_j对应的隐藏状态记为z_j，z_j的取值为区间[1,K]中的一个值，若z_j的取值为1则认为x_j属于第1类隐藏状态，若z_j的取值为k则认为x_j属于第k类隐藏状态，若z_j的取值为K则认为x_j属于第K类隐藏状态；接着计算隐马尔可夫模型中的每个观测数据属于各类隐藏状态的概率，将x_j属于第k类隐藏状态的概率记为

最后计算隐马尔可夫模型中的状态转移概率矩阵，记为Q，

其中，K和k均为正整数，K表示隐马尔可夫模型中的隐藏状态的类别数，2≤K≤10，1≤k≤K，

表示x_j服从的高斯分布的概率密度函数，其变量为x_j、均值为μ_k、方差为τ_k ^-1，μ_k表示属于第k类隐藏状态的高斯分布的均值，τ_k表示属于第k类隐藏状态的高斯分布的精度即方差的倒数，Q_1,1、Q_1,2、Q_1,k'、Q_1,K对应表示Q中的第1行第1列的元素、第1行第2列的元素、第1行第k'列的元素、第1行第K列的元素，Q_2,1、Q_2,2、Q_2,k'、Q_2,K对应表示Q中的第2行第1列的元素、第2行第2列的元素、第2行第k'列的元素、第2行第K列的元素，Q_k,1、Q_k,2、Q_k,k'、Q_k,K对应表示Q中的第k行第1列的元素、第k行第2列的元素、第k行第k'列的元素、第k行第K列的元素，Q_K,1、Q_K,2、Q_K,k'、Q_K,K对应表示Q中的第K行第1列的元素、第K行第2列的元素、第K行第k'列的元素、第K行第K列的元素，1≤k'≤K，Q_k,k'表示z_j'-1＝k的条件下z_j'＝k'的概率，2≤j'≤L，z_j'-1表示隐马尔可夫模型中的第j'-1个观测数据x_j'-1对应的隐藏状态，z_j'表示隐马尔可夫模型中的第j'个观测数据x_j'对应的隐藏状态。

步骤三：在隐马尔可夫模型中引入粘性因子，得到粘性隐马尔可夫模型；在粘性隐马尔可夫模型中，初始化属于每类隐藏状态的高斯分布的均值和精度，将μ_k的初始化值记为

将τ_k的初始化值记为

初始化隐藏状态的类别数K，将K的初始化值记为K⁽⁰⁾，K⁽⁰⁾为区间[2,10]内的任意一个正整数，如K⁽⁰⁾取值为4；初始化状态转移概率矩阵Q，将Q的初始化值记为Q⁽⁰⁾，Q⁽⁰⁾中的每行中的所有元素的共轭先验分布服从狄利克雷分布，Q⁽⁰⁾中的第k⁽⁰⁾行中的所有元素的共轭先验分布服从的狄利克雷分布为：

其中，

表示Q⁽⁰⁾中的第k⁽⁰⁾行中的所有元素，Dir()表示狄利克雷分布，γ表示狄利克雷分布的参数，κ表示粘性因子，δ(k⁽⁰⁾,1)表示两个参数分别为k⁽⁰⁾和1的克罗内克函数，δ(k⁽⁰⁾,k'⁽⁰⁾)表示两个参数分别为k⁽⁰⁾和k'⁽⁰⁾的克罗内克函数，δ(k⁽⁰⁾,K⁽⁰⁾)表示两个参数分别为k⁽⁰⁾和K⁽⁰⁾的克罗内克函数，

γ+κδ(k⁽⁰⁾,1)表示共轭先验分布服从的狄利克雷分布的第1个元素，γ+κδ(k⁽⁰⁾,k'⁽⁰⁾)表示共轭先验分布服从的狄利克雷分布的第k'⁽⁰⁾个元素，γ+κδ(k⁽⁰⁾,K⁽⁰⁾)表示共轭先验分布服从的狄利克雷分布的第K⁽⁰⁾个元素，1≤k⁽⁰⁾≤K⁽⁰⁾，1≤k'⁽⁰⁾≤K⁽⁰⁾。

步骤四：令t表示迭代次数，t的初始值为1；令T表示设定的最大迭代次数，T≥3，在本实施例中取T＝100。

步骤五：计算在第t次迭代下粘性隐马尔可夫模型中的所有观测数据对应的隐藏状态的聚类结果，记为z^(t)，

其中，

表示求使得p(z|x,Q^(t-1),μ^(t-1),τ^(t-1))取最大值时变量z的值，z为粘性隐马尔可夫模型中的所有观测数据对应的隐藏状态构成的向量，z＝[z₁,z₂,…,z_j,…,z_L]，z₁表示第1个观测数据x₁对应的隐藏状态，z₂表示第2个观测数据x₂对应的隐藏状态，z_L表示第L个观测数据x_L对应的隐藏状态，x表示粘性隐马尔可夫模型中的所有观测数据构成的向量，x＝[x₁,x₂,…,x_j,…,x_L]，t＝1时Q^(t-1)即为Q⁽⁰⁾，t≠1时Q^(t-1)表示在第t-1次迭代下粘性隐马尔可夫模型中的状态转移概率矩阵Q的值，t＝1时μ^(t-1)即为μ的初始值μ⁽⁰⁾，μ＝[μ₁,μ₂,…,μ_K]，μ₁表示属于第1类隐藏状态的高斯分布的均值，μ_K表示属于第K类隐藏状态的高斯分布的均值，

表示μ₁的初始化值，

表示μ₂的初始化值，

表示属于第K⁽⁰⁾类隐藏状态的高斯分布的均值

的初始化值，t≠1时μ^(t-1)表示μ在第t-1次迭代下的值，

表示μ₁在第t-1次迭代下的值，

表示μ₂在第t-1次迭代下的值，

表示属于第K^(t-1)类隐藏状态的高斯分布的均值

在第t-1次迭代下的值，t＝1时τ^(t-1)即为τ的初始值τ⁽⁰⁾，τ＝[τ₁,τ₂,…,τ_K]，τ₁表示属于第1类隐藏状态的高斯分布的精度，τ_K表示属于第K类隐藏状态的高斯分布的精度，

表示τ₁的初始化值，

表示τ₂的初始化值，

表示属于第K⁽⁰⁾类隐藏状态的高斯分布的精度

的初始化值，t≠1时τ^(t-1)表示τ在第t-1次迭代下的值，

表示τ₁在第t-1次迭代下的值，

表示τ₂在第t-1次迭代下的值，

表示属于第K^(t-1)类隐藏状态的高斯分布的精度

在第t-1次迭代下的值，t＝1时K^(t-1)即为K的初始化值K⁽⁰⁾，t≠1时K^(t-1)表示K在第t-1次迭代下的值，p(z|x,Q^(t-1),μ^(t-1),τ^(t-1))表示z的后验概率，根据贝叶斯定理得到

p(z_j|x,Q^(t-1),μ^(t-1),τ^(t-1))表示z_j的后验概率，符号“∝”表示正比，x_j+1表示第j+1个观测数据，x_j+2表示第j+2个观测数据，p(z_j,x₁,x₂,...,x_j|Q^(t-1),μ^(t-1),τ^(t-1))表示z_j,x₁,x₂,...,x_j的联合概率，p(x_j+1,x_j+2,...,x_L|z_j,Q^(t-1),μ^(t-1),τ^(t-1))表示z_j的条件下x_j+1,x_j+2,...,x_L的联合概率，p(z_j,x₁,x₂,...,x_j|Q^(t-1),μ^(t-1),τ^(t-1))和p(x_j+1,x_j+2,...,x_L|z_j,Q^(t-1),μ^(t-1),τ^(t-1))通过现有的前向后向算法计算得到。

步骤六：通过生灭过程更新z^(t)，进而计算在第t次迭代下粘性隐马尔可夫模型中的隐藏状态的类别数K的值K^(t)，具体过程为：

1)统计z^(t)中值等于区间[1,K^(t-1)]中的每个值的元素的总个数，将z^(t)中值等于k^(t-1)的元素的总个数记为

2)按1,…,k^(t-1),…,K^(t-1)从小到大的顺序排列Num₁至

得到统计个数排列序列。

3)若统计个数排列序列中只有一个0值且K^(t-1)≠2，则假设该0值对应区间[1,K^{(t -1)}]中的k^(t-1)，那么将z^(t)中值分别等于(k+1)^(t-1)至K^(t-1)的所有元素的值均减1，将更新后的z^(t)重新记为z^*(t)，并令K^(t)＝K^(t-1)-1，再执行步骤七；若统计个数排列序列中只有一个0值且K^(t-1)＝2，或者统计个数排列序列中有多个0值，则执行步骤4)；若统计个数排列序列中没有0值且K^(t-1)＜10，则从1到L中随机产生一个起始值，记为L_min，从L_min到L中随机产生一个终止值，记为L_max，将z^(t)中的

的值均设为K^(t-1)+1，将更新后的z^(t)重新记为z^*(t)，并令K^(t)＝K^(t-1)+1，再执行步骤七；如果统计个数排列序列为除上述四种情况外的其余情况，则保持z^(t)不变，将z^(t)重新记为z^*(t)，并令K^(t)＝K^(t-1)，再执行步骤七。

4)当统计个数排列序列中的第ω个值为L而第1个至第ω-1个以及第ω+1个至第K^(t-1)个值均为0即z^(t)中的所有元素的值等于ω时，从1到L中随机产生一个起始值，记为L'_min，从L'_min到L中随机产生一个终止值，记为L'_max，将z^(t)中的

的值均设为2，将z^(t)中的

和

的值均设为1，将更新后的z^(t)重新记为z^*(t)，并令K^(t)＝2，再执行步骤七；当统计个数排列序列中有多个0值且任一个非0值不为L时，执行步骤5)。

5)设0值有ξ个，针对第1个0值和第2个0值，假设第1个0值对应区间[1,K^(t-1)]中的k^(t-1)、第2个0值对应区间[1,K^(t-1)]中的(k+υ)^(t-1)，那么将z^(t)中值分别等于(k+1)^(t-1)至(k+υ-1)^(t-1)的所有元素的值均减1；针对第2个0值和第3个0值，假设第2个0值对应区间[1,K^{(t -1)}]中的(k+υ)^(t-1)、第3个0值对应区间[1,K^(t-1)]中的

那么将z^(t)中值分别等于(k+υ+1)^(t-1)至

的所有元素的值均减2；依次类推，针对第ξ-1个0值和第ξ个0值，假设第ξ-1个0值对应区间[1,K^(t-1)]中的(k+θ)^(t-1)、第ξ个0值对应区间[1,K^(t-1)]中的(k+ρ)^(t-1)，那么将z^(t)中值分别等于(k+θ+1)^(t-1)至(k+ρ-1)^(t-1)的所有元素的值均减ξ-1；而后将z^(t)中值分别等于(k+ρ+1)^(t-1)至K^(t-1)的所有元素的值均减ξ；将更新后的z^(t)重新记为z^*(t)，并令K^(t)＝K^(t-1)-ξ，再执行步骤七。

上述，k^(t-1)为区间[1,K^(t-1)]中的第k^(t-1)个值，Num₁表示z^(t)中值等于1的元素的总个数，

表示z^(t)中值等于K^(t-1)的元素的总个数，1≤ω≤K^(t-1)，1≤L_min≤L_max≤L，

对应表示z^(t)中的第L_min个元素、第L_min+1个元素、…、第L_max个元素，1≤L'_min≤L'_max≤L，

对应表示z^(t)中的第L'_min个元素、第L'_min+1个元素、…、第L'_max个元素，

对应表示z^(t)中的第1个元素、第2个元素、…、第L'_min-1个元素，

对应表示z^(t)中的第L'_max+1个元素、第L'_max+2个元素、…、第L个元素，1＜ξ＜K^(t-1)，

K^(t)＝K^(t-1)中的“＝”为赋值符号，

表示z₁在第t次迭代下经过生灭过程后的值，

表示z₂在第t次迭代下经过生灭过程后的值，

表示z_j在第t次迭代下经过生灭过程后的值，

表示z_L在第t次迭代下经过生灭过程后的值。

该过程又更新了z^(t)，该过程为生灭过程。

步骤七：计算在第t次迭代下粘性隐马尔可夫模型中的状态转移概率矩阵Q的值，记为Q^(t)，Q^(t)中的第k^(t)行中的所有元素的共轭先验分布服从的狄利克雷分布为：

Q^(t)中的第k^(t)行中的所有元素的后验分布服从的狄利克雷分布为：

其中，1≤k^(t)≤K^(t)，1≤k'^(t)≤K^(t)，

表示Q^(t)中的第k^(t)行中的所有元素，

表示在第t次迭代下经过生灭过程后从第k^(t)类隐藏状态转移到第1类隐藏状态的观测数据的数量，

表示在第t次迭代下经过生灭过程后从第k^(t)类隐藏状态转移到第k'^(t)类隐藏状态的观测数据的数量，

表示在第t次迭代下经过生灭过程后从第k^(t)类隐藏状态转移到第K^(t)类隐藏状态的观测数据的数量，δ(k^(t),1)表示两个参数分别为k^(t)和1的克罗内克函数，δ(k^(t),k'^(t))表示两个参数分别为k^(t)和k'^(t)的克罗内克函数，δ(k^(t),K^(t))表示两个参数分别为k^(t)和K^(t)的克罗内克函数，

表示后验分布服从的狄利克雷分布的第1个元素，

表示后验分布服从的狄利克雷分布的第k'^(t)个元素，

表示后验分布服从的狄利克雷分布的第K^(t)个元素。

步骤八：利用属于同一类隐藏状态的所有观测数据，根据贝叶斯定理，计算在x和z^*(t)确定后μ在第t次迭代下的值μ^(t)和τ在第t次迭代下的值τ^(t)的后验概率，记为p(μ^(t),τ^(t)|x,z^*(t))，

其中，μ^(t)表示μ在第t次迭代下的值，

表示μ₁在第t次迭代下的值，

表示属于第k^(t)类隐藏状态的高斯分布的均值

在第t次迭代下的值，

表示属于第K^(t)类隐藏状态的高斯分布的均值

在第t次迭代下的值，τ^(t)表示τ在第t次迭代下的值，

表示τ₁在第t次迭代下的值，

表示属于第k^(t)类隐藏状态的高斯分布的精度

在第t次迭代下的值，

表示属于第K^(t)类隐藏状态的高斯分布的精度

在第t次迭代下的值，

表示在第t次迭代下经过生灭过程后属于第k^(t)类隐藏状态的观测数据的数量，

表示在第t次迭代下经过生灭过程后属于第k^(t)类隐藏状态的所有观测数据的平均值，

表示在第t次迭代下经过生灭过程后属于第k^(t)类隐藏状态的第

个观测数据，

表示

服从的高斯分布的概率密度函数，其变量为

均值为

方差为

表示

服从的伽马分布的概率密度函数，其变量为

形状参数为

速率参数为

η₀、m₀、a₀、b₀均为常数，在本实施例中取η₀＝1、m₀＝1、a₀＝1、b₀＝1。

步骤九：判断t＜T是否成立，如果成立，则令t＝t+1，然后返回步骤五继续迭代；如果不成立，则执行步骤十；其中，t＝t+1中的“＝”为赋值符号。

步骤十：将μ^(t)中的每个元素的值作为对应一类隐藏状态的功率估计值，即将μ^(t)中的第k^(t)个元素的值作为第k^(tw类隐藏状态的功率估计值；然后根据每类隐藏状态的功率估计值和每个观测数据对应的隐藏状态在第t次迭代下经过生灭过程后的值，计算每个频带的功率估计值，将第j个频带的功率估计值记为β_j，若

则β_j等于第k^(t)类隐藏状态的功率估计值；再对每个频带的功率估计值与门限值进行比较，对于β_j，如果β_j小于门限值，则认为第j个频带未被授权用户占用，并将第j个频带作为可用频带；如果β_j大于或等于门限值，则认为第j个频带已被授权用户占用，不能被认知用户使用；其中，门限值是根据给定的虚警概率计算得到的，在本实施例中给定的虚警概率为0.1。在具体实施时也可直接将门限值设定为4.0548，该值是通过大量实验获得的。

通过以下仿真来进一步说明本发明方法的可行性和有效性。

在仿真中，噪声功率为

最大迭代次数T＝100，m₀、η₀、a₀、b₀均取1，狄利克雷分布的参数γ＝1，粘性因子κ＝50，样本数量N＝1000，蒙特卡洛次数为1000。

图2为信噪比为-10dB，频带占用率为20％(频带具有稀疏性)时本发明方法和使用无粘性的隐马尔可夫模型的方法的平均功率估计误差随迭代次数的变化曲线。平均功率估计误差的计算公式为：

其中β_j表示第j个频带的功率估计值，

表示第j个频带内授权用户的信号功率，

表示噪声功率，符号

表示求2-范数的平方符号。从图2中可以看出，本发明方法收敛速度快，并且与使用无粘性的隐马尔可夫模型的方法相比，具有更低的平均功率估计误差。

图3为信噪比为-10dB，频带占用率为20％(频带具有稀疏性)时本发明方法和使用无粘性的隐马尔可夫模型的方法的ROC曲线。图4为信噪比为-10dB，频带占用率为80％(频带不具有稀疏性)时本发明方法和使用无粘性的隐马尔可夫模型的方法的ROC曲线。从图3和图4中可以看出，本发明方法不管频带是否具有稀疏性，均能够对多个频带进行频谱感知，并且与使用无粘性的隐马尔可夫模型的方法相比，具有良好的检测性能。因此本发明方法确实是能够很好地提高多频带频谱感知的检测概率。

Claims (1)

1.一种基于生灭过程和粘性隐马尔可夫模型的多带频谱感知方法，其特征在于具体包括以下步骤：

步骤一：在认知无线电系统中，对连续的L个频带内的信号进行采样，且对每个频带内的信号进行等时间间隔地采样，共采样得到N个样本，将对第j个频带内的信号进行采样得到的第n个样本记为r_j(n)；然后计算每个频带对应的接收信号功率，将第j个频带对应的接收信号功率记为x_j，即为对第j个频带内的信号进行采样得到的所有样本的平均功率，

其中，L、N、j和n均为正整数，L＞1，100≤N≤1000，1≤j≤L，1≤n≤N，符号“||”为求模符号，x_j服从高斯分布，即

表示噪声功率，

表示第j个频带内授权用户的信号功率，

表示第j个频带未被授权用户占用，

表示第j个频带已被授权用户占用，

表示x_j服从均值为

方差为

的高斯分布，

表示x_j服从均值为

方差为

的高斯分布；

步骤二：将每个频带对应的接收信号功率作为隐马尔可夫模型中的观测数据，即隐马尔可夫模型中的第j个观测数据为x_j；然后确定隐马尔可夫模型中的每个观测数据所对应的一个隐藏状态，将x_j对应的隐藏状态记为z_j，z_j的取值为区间[1,K]中的一个值，若z_j的取值为1则认为x_j属于第1类隐藏状态，若z_j的取值为k则认为x_j属于第k类隐藏状态，若z_j的取值为K则认为x_j属于第K类隐藏状态；接着计算隐马尔可夫模型中的每个观测数据属于各类隐藏状态的概率，将x_j属于第k类隐藏状态的概率记为

最后计算隐马尔可夫模型中的状态转移概率矩阵，记为Q，

其中，K和k均为正整数，K表示隐马尔可夫模型中的隐藏状态的类别数，2≤K≤10，1≤k≤K，

表示x_j服从的高斯分布的概率密度函数，其变量为x_j、均值为μ_k、方差为

μ_k表示属于第k类隐藏状态的高斯分布的均值，τ_k表示属于第k类隐藏状态的高斯分布的精度即方差的倒数，Q_1,1、Q_1,2、Q_1,k'、Q_1,K对应表示Q中的第1行第1列的元素、第1行第2列的元素、第1行第k'列的元素、第1行第K列的元素，Q_2,1、Q_2,2、Q_2,k'、Q_2,K对应表示Q中的第2行第1列的元素、第2行第2列的元素、第2行第k'列的元素、第2行第K列的元素，Q_k,1、Q_k,2、Q_k,k'、Q_k,K对应表示Q中的第k行第1列的元素、第k行第2列的元素、第k行第k'列的元素、第k行第K列的元素，Q_K,1、Q_K,2、Q_K,k'、Q_K,K对应表示Q中的第K行第1列的元素、第K行第2列的元素、第K行第k'列的元素、第K行第K列的元素，1≤k'≤K，Q_k,k'表示z_j'-1＝k的条件下z_j'＝k'的概率，2≤j'≤L，z_j'-1表示隐马尔可夫模型中的第j'-1个观测数据x_j'-1对应的隐藏状态，z_j'表示隐马尔可夫模型中的第j'个观测数据x_j'对应的隐藏状态；

步骤三：在隐马尔可夫模型中引入粘性因子，得到粘性隐马尔可夫模型；在粘性隐马尔可夫模型中，初始化属于每类隐藏状态的高斯分布的均值和精度，将μ_k的初始化值记为

将τ_k的初始化值记为

初始化隐藏状态的类别数K，将K的初始化值记为K⁽⁰⁾，K⁽⁰⁾为区间[2,10]内的任意一个正整数；初始化状态转移概率矩阵Q，将Q的初始化值记为Q⁽⁰⁾，Q⁽⁰⁾中的每行中的所有元素的共轭先验分布服从狄利克雷分布，Q⁽⁰⁾中的第k⁽⁰⁾行中的所有元素的共轭先验分布服从的狄利克雷分布为：

其中，

表示Q⁽⁰⁾中的第k⁽⁰⁾行中的所有元素，Dir()表示狄利克雷分布，γ表示狄利克雷分布的参数，κ表示粘性因子，δ(k⁽⁰⁾,1)表示两个参数分别为k⁽⁰⁾和1的克罗内克函数，δ(k⁽⁰⁾,k'⁽⁰⁾)表示两个参数分别为k⁽⁰⁾和k'⁽⁰⁾的克罗内克函数，δ(k⁽⁰⁾,K⁽⁰⁾)表示两个参数分别为k⁽⁰⁾和K⁽⁰⁾的克罗内克函数，

γ+κδ(k⁽⁰⁾,1)表示共轭先验分布服从的狄利克雷分布的第1个元素，γ+κδ(k⁽⁰⁾,k'⁽⁰⁾)表示共轭先验分布服从的狄利克雷分布的第k'⁽⁰⁾个元素，γ+κδ(k⁽⁰⁾,K⁽⁰⁾)表示共轭先验分布服从的狄利克雷分布的第K⁽⁰⁾个元素，1≤k⁽⁰⁾≤K⁽⁰⁾，1≤k'⁽⁰⁾≤K⁽⁰⁾；

步骤四：令t表示迭代次数，t的初始值为1；令T表示设定的最大迭代次数，T≥3；

步骤五：计算在第t次迭代下粘性隐马尔可夫模型中的所有观测数据对应的隐藏状态的聚类结果，记为z^(t)，

其中，

表示求使得p(z|x,Q^(t-1),μ^(t-1),τ^(t-1))取最大值时变量z的值，z为粘性隐马尔可夫模型中的所有观测数据对应的隐藏状态构成的向量，z＝[z₁,z₂,…,z_j,…,z_L]，z₁表示第1个观测数据x₁对应的隐藏状态，z₂表示第2个观测数据x₂对应的隐藏状态，z_L表示第L个观测数据x_L对应的隐藏状态，x表示粘性隐马尔可夫模型中的所有观测数据构成的向量，x＝[x₁,x₂,…,x_j,…,x_L]，t＝1时Q^(t-1)即为Q⁽⁰⁾，t≠1时Q^(t-1)表示在第t-1次迭代下粘性隐马尔可夫模型中的状态转移概率矩阵Q的值，t＝1时μ^(t-1)即为μ的初始值μ⁽⁰⁾，μ＝[μ₁,μ₂,…,μ_K]，μ₁表示属于第1类隐藏状态的高斯分布的均值，μ_K表示属于第K类隐藏状态的高斯分布的均值，

表示μ₁的初始化值，

表示μ₂的初始化值，

表示属于第K⁽⁰⁾类隐藏状态的高斯分布的均值

的初始化值，t≠1时μ^(t-1)表示μ在第t-1次迭代下的值，

表示μ₁在第t-1次迭代下的值，

表示μ₂在第t-1次迭代下的值，

表示属于第K^(t-1)类隐藏状态的高斯分布的均值

在第t-1次迭代下的值，t＝1时τ^(t-1)即为τ的初始值τ⁽⁰⁾，τ＝[τ₁,τ₂,…,τ_K]，τ₁表示属于第1类隐藏状态的高斯分布的精度，τ_K表示属于第K类隐藏状态的高斯分布的精度，

表示τ₁的初始化值，

表示τ₂的初始化值，

表示属于第K⁽⁰⁾类隐藏状态的高斯分布的精度

的初始化值，t≠1时τ^(t-1)表示τ在第t-1次迭代下的值，

表示τ₁在第t-1次迭代下的值，

表示τ₂在第t-1次迭代下的值，

表示属于第K^(t-1)类隐藏状态的高斯分布的精度

在第t-1次迭代下的值，t＝1时K^(t-1)即为K的初始化值K⁽⁰⁾，t≠1时K^(t-1)表示K在第t-1次迭代下的值，p(z|x,Q^(t-1),μ^(t-1),τ^(t-1))表示z的后验概率，根据贝叶斯定理得到

p(z_j|x,Q^(t-1),μ^(t-1),τ^(t-1))表示z_j的后验概率，符号“∝”表示正比，x_j+1表示第j+1个观测数据，x_j+2表示第j+2个观测数据，p(z_j,x₁,x₂,...,x_j|Q^(t-1),μ^(t-1),τ^(t-1))表示z_j,x₁,x₂,...,x_j的联合概率，p(x_j+1,x_j+2,...,x_L|z_j,Q^(t-1),μ^(t-1),τ^(t-1))表示z_j的条件下x_j+1,x_j+2,...,x_L的联合概率，p(z_j,x₁,x₂,...,x_j|Q^(t-1),μ^(t-1),τ^(t-1))和p(x_j+1,x_j+2,...,x_L|z_j,Q^(t-1),μ^(t-1),τ^(t-1))通过前向后向算法计算得到；

步骤六：通过生灭过程更新z^(t)，进而计算在第t次迭代下粘性隐马尔可夫模型中的隐藏状态的类别数K的值K^(t)，具体过程为：

1)统计z^(t)中值等于区间[1,K^(t-1)]中的每个值的元素的总个数，将z^(t)中值等于k^(t-1)的元素的总个数记为

2)按1,…,k^(t-1),…,K^(t-1)从小到大的顺序排列Num₁至

得到统计个数排列序列；

3)若统计个数排列序列中只有一个0值且K^(t-1)≠2，则假设该0值对应区间[1,K^(t-1)]中的k^(t-1)，那么将z^(t)中值分别等于(k+1)^(t-1)至K^(t-1)的所有元素的值均减1，将更新后的z^(t)重新记为z^*(t)，并令K^(t)＝K^(t-1)-1，再执行步骤七；若统计个数排列序列中只有一个0值且K^(t-1)＝2，或者统计个数排列序列中有多个0值，则执行步骤4)；若统计个数排列序列中没有0值且K^(t-1)＜10，则从1到L中随机产生一个起始值，记为L_min，从L_min到L中随机产生一个终止值，记为L_max，将z^(t)中的

的值均设为K^(t-1)+1，将更新后的z^(t)重新记为z^*(t)，并令K^(t)＝K^(t-1)+1，再执行步骤七；如果统计个数排列序列为除上述四种情况外的其余情况，则保持z^(t)不变，将z^(t)重新记为z^*(t)，并令K^(t)＝K^(t-1)，再执行步骤七；

4)当统计个数排列序列中的第ω个值为L而第1个至第ω-1个以及第ω+1个至第K^(t-1)个值均为0即z^(t)中的所有元素的值等于ω时，从1到L中随机产生一个起始值，记为L'_min，从L'_min到L中随机产生一个终止值，记为L'_max，将z^(t)中的

的值均设为2，将z^(t)中的

和

的值均设为1，将更新后的z^(t)重新记为z^*(t)，并令K^(t)＝2，再执行步骤七；当统计个数排列序列中有多个0值且任一个非0值不为L时，执行步骤5)；

5)设0值有ξ个，针对第1个0值和第2个0值，假设第1个0值对应区间[1,K^(t-1)]中的k^(t-1)、第2个0值对应区间[1,K^(t-1)]中的(k+υ)^(t-1)，那么将z^(t)中值分别等于(k+1)^(t-1)至(k+υ-1)^(t-1)的所有元素的值均减1；针对第2个0值和第3个0值，假设第2个0值对应区间[1,K^(t-1)]中的(k+υ)^(t-1)、第3个0值对应区间[1,K^(t-1)]中的

那么将z^(t)中值分别等于(k+υ+1)^(t-1)至

的所有元素的值均减2；依次类推，针对第ξ-1个0值和第ξ个0值，假设第ξ-1个0值对应区间[1,K^(t-1)]中的(k+θ)^(t-1)、第ξ个0值对应区间[1,K^(t-1)]中的(k+ρ)^(t-1)，那么将z^(t)中值分别等于(k+θ+1)^(t-1)至(k+ρ-1)^(t-1)的所有元素的值均减ξ-1；而后将z^(t)中值分别等于(k+ρ+1)^(t-1)至K^(t-1)的所有元素的值均减ξ；将更新后的z^(t)重新记为z^*(t)，并令K^(t)＝K^(t-1)-ξ，再执行步骤七；

上述，k^(t-1)为区间[1,K^(t-1)]中的第k^(t-1)个值，Num₁表示z^(t)中值等于1的元素的总个数，

表示z^(t)中值等于K^(t-1)的元素的总个数，1≤ω≤K^(t-1)，1≤L_min≤L_max≤L，

对应表示z^(t)中的第L_min个元素、第L_min+1个元素、…、第L_max个元素，1≤L'_min≤L'_max≤L，

对应表示z^(t)中的第L'_min个元素、第L'_min+1个元素、…、第L'_max个元素，

对应表示z^(t)中的第1个元素、第2个元素、…、第L'_min-1个元素，

对应表示z^(t)中的第L'_max+1个元素、第L'_max+2个元素、…、第L个元素，1＜ξ＜K^(t-1)，

K^(t)＝K^(t-1)中的“＝”为赋值符号，

表示z₁在第t次迭代下经过生灭过程后的值，

表示z₂在第t次迭代下经过生灭过程后的值，

表示z_j在第t次迭代下经过生灭过程后的值，

表示z_L在第t次迭代下经过生灭过程后的值；

步骤七：计算在第t次迭代下粘性隐马尔可夫模型中的状态转移概率矩阵Q的值，记为Q^(t)，Q^(t)中的第k^(t)行中的所有元素的共轭先验分布服从的狄利克雷分布为：

Q^(t)中的第k^(t)行中的所有元素的后验分布服从的狄利克雷分布为：

；其中，1≤k^(t)≤K^(t)，1≤k'^(t)≤K^(t)，

表示Q^(t)中的第k^(t)行中的所有元素，

表示在第t次迭代下经过生灭过程后从第k^(t)类隐藏状态转移到第1类隐藏状态的观测数据的数量，

表示在第t次迭代下经过生灭过程后从第k^(t)类隐藏状态转移到第k'^(t)类隐藏状态的观测数据的数量，

表示在第t次迭代下经过生灭过程后从第k^(t)类隐藏状态转移到第K^(t)类隐藏状态的观测数据的数量，δ(k^(t),1)表示两个参数分别为k^(t)和1的克罗内克函数，δ(k^(t),k'^(t))表示两个参数分别为k^(t)和k'^(t)的克罗内克函数，δ(k^(t),K^(t))表示两个参数分别为k^(t)和K^(t)的克罗内克函数，

表示后验分布服从的狄利克雷分布的第1个元素，

表示后验分布服从的狄利克雷分布的第k'^(t)个元素，

表示后验分布服从的狄利克雷分布的第K^(t)个元素；

步骤八：利用属于同一类隐藏状态的所有观测数据，根据贝叶斯定理，计算在x和z^*(t)确定后μ在第t次迭代下的值μ^(t)和τ在第t次迭代下的值τ^(t)的后验概率，记为p(μ^(t),τ^(t)|x,z^*(t))，

其中，μ^(t)表示μ在第t次迭代下的值，

表示μ₁在第t次迭代下的值，

表示属于第k^(t)类隐藏状态的高斯分布的均值

在第t次迭代下的值，

表示属于第K^(t)类隐藏状态的高斯分布的均值

在第t次迭代下的值，τ^(t)表示τ在第t次迭代下的值，

表示τ₁在第t次迭代下的值，

表示属于第k^(t)类隐藏状态的高斯分布的精度

在第t次迭代下的值，

表示属于第K^(t)类隐藏状态的高斯分布的精度

在第t次迭代下的值，

表示在第t次迭代下经过生灭过程后属于第k^(t)类隐藏状态的观测数据的数量，

表示在第t次迭代下经过生灭过程后属于第k^(t)类隐藏状态的所有观测数据的平均值，

表示在第t次迭代下经过生灭过程后属于第k^(t)类隐藏状态的第

个观测数据，

表示

服从的高斯分布的概率密度函数，其变量为

均值为

方差为

表示

服从的伽马分布的概率密度函数，其变量为

形状参数为

速率参数为

η₀、m₀、a₀、b₀均为常数；

步骤九：判断t＜T是否成立，如果成立，则令t＝t+1，然后返回步骤五继续迭代；如果不成立，则执行步骤十；其中，t＝t+1中的“＝”为赋值符号；

步骤十：将μ^(t)中的每个元素的值作为对应一类隐藏状态的功率估计值，即将μ^(t)中的第k^(t)个元素的值作为第k^(t)类隐藏状态的功率估计值；然后根据每类隐藏状态的功率估计值和每个观测数据对应的隐藏状态在第t次迭代下经过生灭过程后的值，计算每个频带的功率估计值，将第j个频带的功率估计值记为β_j，若

则β_j等于第k^(t)类隐藏状态的功率估计值；再对每个频带的功率估计值与门限值进行比较，对于β_j，如果β_j小于门限值，则认为第j个频带未被授权用户占用，并将第j个频带作为可用频带；如果β_j大于或等于门限值，则认为第j个频带已被授权用户占用，不能被认知用户使用；其中，门限值是根据给定的虚警概率计算得到的。

Images