基于贝叶斯信息准则的信源数估计方法.pdf

本发明提供了一种在贝叶斯信息准则(BIC)框架下的信源数估计方法，适用于大规模自适应天线场景，能够在广义渐近条件下，即m，n→∞，m/n→c∈(0，∞)，m和n分别代表天线数和快拍数，在该条件下提供可靠的信源数目检测。本发明由对数似然函数和代价函数共同计算获得先验概率，并通过最大化先验概率，有效地得到了信源的数目。仿真结果证明本发明的信源数估计方法的优越性以及有效性。

1.  一种大规模自适应天线下基于贝叶斯信息准则的信源数估计方法，适用 于m阵元的均匀线阵，有d个远场窄带信号{s₁(t)，L，s_d(t)}从不同方向 入射到所述阵列，假设信源和阵列处在同一平面内，则第t个 快拍时，所述阵列输出可以表示为x_t＝As_t+w_t，(t＝1，L，n)，其中， x_t＝[x₁(t)，L，x_mt()^T]∈？^m×1，s_t＝[s₁(t)，L，s_d(t)]^T∈？^d×1， w_t＝[w₁(t)，L，w_m(t)]^T∈？^m×1分别表示观测向量、阵列流形、信号矢量和噪 声矢量，d是第i个信源的波达方向对应的导向矢量，(·)^T表 示转置运算，d是未知信源数目，m是天线数目，n是快拍数；其特征 在于：所述方法包括以下步骤：
步骤1：计算采样协方差矩阵令1₁，L，1_m和e₁，L，e_m，分别为递减特征值和对应的特征向量；
步骤2：执行的特征值分解，得到1₁，L，1_m；
步骤3：利用1₁，L，1_m计算a(k)、g(k)、和P(k，m，n)，其中：
τ ^ k = ( 1 / ( m - k ) ) Σ i = k + 1 m 1 i , a ( k ) = 1 m - k Σ i = k + 1 m 1 i , g ( k ) = ( Π i = k + 1 m 1 i ) 1 m - k ]]>分别是算术平均值 和几何平均值， P ( k , m , n ) @ mk ( log ( 2 n ) - 1 k Σ i = 1 k log 1 i τ ^ k ) , ]]>0≤k≤m-1；
步骤4：计算m个贝叶斯值，0≤k≤m-1：
BIC ( k ) = 2 n ( m - k ) log a ( k ) g ( k ) + P ( k , m , n ) ]]>
步骤5：最小化准则下式得到信源数估计，其中，为信源数估计值：
d ^ = arg min k = 0 , L , m ‾ - 1 BIC ( k ) . ]]>

基于贝叶斯信息准则的信源数估计方法
技术领域
本发明涉及信号处理领域，特别涉及一种大规模自适应天线下的信源 数估计方法。
背景技术
对频谱利用率的突出贡献，使得大规模MIMO系统受到了广泛的关 注。大规模MIMO系统的基站(BS)使用了大量的天线来进行发射和接收， 这使传统基于子空间的波达方向(DOA)估计算法遭受了严重的性能下降。 在天线数量与采样数量相当时，基于子空间的算法将无法求得正确的子空 间，为了解决这个问题，目前提出了许多更有效的子空间算法来处理大规 模MIMO问题，然而，信源数估计作为子空间算法的前提，目前却没有任 何相应的算法，在低信噪比(SNR)和小样本环境下这是一个巨大的挑战。
传统的信源数估计方法中，假设检验包含的球形检验和随机矩阵理论 检验(RMT)都需要寻找主观阈值。其中RMT方法能够提供比经典最小描述 长度算法(MDL)更小的检测阈值。区别于假设检验方法，信息理论准则法 (ITC)，如赤池信息量准则(AIC)，贝叶斯信息准则(BIC)法是从信息论角度 推导而出，并且不需要用户自定义参数。为处理相干信号，结合了最大似 然估计(ML)和MDL准则的联合方法，作为增强的MDL准则能有效处理 相干源信号。另外，优于MDL方法的预测描述长度(PDL)算法也可以处理 相干信号，但每个快拍需要计算最大似然函数而导致计算复杂度较高。在 非均匀噪声环境下，一种利用多级维纳滤波器的最小均方误差算法解决了 非均匀噪声对采样特征值的干扰，改善了计算复杂和鲁棒性信源数估计的 问题。
上述算法都是假设天线数m固定，且快拍数趋于无穷的经典渐进条件 下提出的。实际上，广义渐近条件，即m，n→∞，m/n→c∈(0，∞)，更适用 于天线数与快拍数都较大的大规模MIMO应用。另一方面，广义渐近条件 能够对快拍数和天线数有限且数量级相当的实际情况提供更精确的描述。 目前，围绕广义渐近条件的信源数估计，DOA估计和波束形成问题已展开 了大量的科研工作。
ITC方法需要求解KL信息最小化问题，但该问题要求天线数固定且 快拍数趋于无穷，这意味着ITC不适用于广义渐近条件，为改进该算法， 提出的RMT-AIC准则能够在广义渐近条件下正确检测信源数目，但是不 能保证信源数估计为一致估计。为了改善先行回归模型只能处理小样本情 况，改进的AIC方法通过对KL信息的Bootstrap估计进行渐进近似得到。 由此看出，大规模MIMO中的信源数估计是十分重要的，尤其在广义渐近 条件下探索信源数估计的一致方法是相当有意义的。
人们更想得到能在大快拍数情况下估计真实信源数目的信源数估计方 法。相比AIC方法，BIC方法提供了强一致性的保证，传统的BIC准则由 最大似然函数(LF)和代价函数(PF)组成，分别对应数据拟合和模型复杂度。 BIC准则的最小化实际上是数据拟合和模型复杂度的权衡过程，从而得到 模型阶次或信源数目估计。如上所述，已有的BIC准则能最小化一般模型 和拟合逼近模型之间的相关KL信息，却只适用于m固定，n→∞的情况。 在广义渐近条件下，经典BIC准则不能保证估计信源数目的正确性。
发明内容
为了解决上述问题，本发明提出了适用于广义渐近条件的BIC准则。 特别地，通过渐进方法计算广义渐近条件下的LF和PF，得到改进的BIC 准则，并估计信源数目。
本发明采取了以下技术方案：
一种大规模自适应天线下基于贝叶斯信息准则的信源数估计方法，适 用于m阵元的均匀线阵，有d个远场窄带信号(s₁(t)，L，s_d(t)}从不同方向 入射到所述阵列，假设信源和阵列处在同一平面内，则第t个快拍 时，所述阵列输出可以表示为x_t＝As_t+w_t，(t＝1，L，n)，其中， x_t＝[x₁(t)，L，x_m(t)]^T∈？^m×1，s_t＝[s₁(t)，L，s_d(t)]^T∈？^d×1， w_t＝[w₁(t)，L，w_m(t)]^T∈？^m×1分别表示观测向量、阵列流形、信号矢量和噪声矢 量，d是第i个信源的波达方向对应的导向矢量，(·)^T表示转置 运算，d是未知信源数目，m是天线数目，n是快拍数；其特征在于：所述 方法包括以下步骤：
步骤1：计算采样协方差矩阵令l₁，L，l_m和e₁，L，e_m，分别为递减特征值和对应的特征向量；
步骤2：执行的特征值分解，得到l₁，L，l_m；
步骤3：利用l₁，L，l_m计算a(k)、g(k)、和P(k，m，n)，其中：
τ ^ k = ( 1 / ( m - k ) ) Σ i = k + 1 m l i , a ( k ) = 1 m - k Σ i = k + 1 m l i , g ( k ) = ( Π i = k + 1 m l i ) 1 m - k ]]>分别是算术平均值 和几何平均值，
步骤4：计算m个贝叶斯值，0≤k≤m-1：
BIC ( k ) = 2 n ( m - k ) log a ( k ) g ( k ) + P ( k , m , n ) ]]>
步骤5：最小化准则下式得到信源数估计，其中，为信源数估计值：
d ^ = arg min k = 0 , L , m ‾ - 1 BIC ( k ) . ]]>
本发明的有益效果是：本发明提出的方法适用于广义渐近条件的BIC 准则。特别地，通过渐进方法计算广义渐近条件下的LF和PF，得到改进 的BIC准则，并估计信源数目。本发明的方法为大规模MIMO系统中DOA 估计和波束形成确定信号和噪声的子空间提供了精确的保证。
附图说明
图1是本发明的基于贝叶斯信息准则的信源数估计方法的流程图；
图2(a)是本发明的方法和现有技术的算法的正确检测概率与信噪比 的关系(n＝60)示意图；
图2(b)是本发明的方法和现有技术的算法的正确检测概率与信噪比 的关系(n＝150)示意图；
图3(a)是本发明的方法和现有技术的算法的正确检测概率与DOA 间距的关系(n＝60)示意图；
图3(b)是本发明的方法和现有技术的算法的正确检测概率与DOA 间距的关系(n＝150)示意图；
图4(a)是本发明的方法和现有技术的算法的正确检测概率与信噪比 的关系(d＝3)示意图；
图4(b)是本发明的方法和现有技术的算法的正确检测概率与信噪比 的关系(d＝8)示意图；
图5(a)是本发明的方法和现有技术的算法的正确检测概率与快拍数 的关系(d＝3)示意图；
图5(b)是本发明的方法和现有技术的算法的正确检测概率与快拍数 的关系(d＝8)示意图。
具体实施方式
下面结合附图说明及具体实施方式对本发明进一步说明。
考虑一m阵元的均匀线阵，有d个远场窄带信号(s₁(t)，L，s_d(t)}从不同方 向入射到该阵列，假设信源和阵列处在同一平面内，则第t个快拍 时，阵列输出可以表示为
x_t＝As_t+w_t，(t＝1，L，n)    (1)
式中，x_t＝[x₁(t)，L，x_m(t)]^T∈？^m×1，s_t＝[s₁(t)，L，s_d(t)]^T∈？^d×1， w_t＝[w₁(t)，L，w_m(t)]^T∈？^m×1分别表示观测向量、阵列流形(steering matrix)、 信号矢量和噪声矢量。其中，d是第i个信源的波达方向对应 的导向矢量，(·)^T表示转置运算，d是未知信源数目，m是天线数目，n是 快拍数。为简单起见，假设m＜n，信源数固定且小于常数远小于 min(m，n)，即非相干信号独立同分布，即s_t～CN(0_d，R_s)，其 中0_d是d×1零向量，是满秩矩阵，(·)^H表示共轭转置，E[·]是 数学期望。CN(ν，R)表示均值为ν，协方差为R的复高斯分布，～表示“服 从于”。假设噪声为零均值，协方差为τI_m的IID复高斯向量，即 n_t～CN(0_m，τI_m)，其中I_m是m×m单位矩阵，与信号相互独立。
基于上述假设，观测样本可以看作IID高斯向量，x_t～CN(o_m，R)，R是 总体协方差矩阵
R = E [ x t x t H ] = A R s A H + τ I m . - - - ( 2 ) ]]>
由于信号非相干且d＜m，所以R_s非奇异，A列满秩。不失一般性，假设R 的总体特征向量表示为
λ₁≥L≥λ_d≥λ_d+1＝L＝λ_m＝τ.    (3)
相应的总体特征向量表示为u₁，L，u_m，给定公式(3)可直接利用乘子τ去确定 信号的个数。实际上，只能获得采样协方差矩阵

令l₁，L，l_m和e₁，L，e_m，分别为递减特征值和对应的特征向量。因此，在 m，n→∞和m/n→c环境下，本发明的主要工作是从嘈杂的观测(x₁，L，x_n}中 推断出信源数目d。
IID高斯观测信号X＝[x₁，L，x_n]，其联合概率密度函数(PDF)为
f ( X | θ ) = Π t = 1 n 1 π m | R | exp ( - x t H R - 1 x t ) - - - ( 5 ) ]]>
其中|·|表示行列式，θ是真实模型的未知向量参数，给定为 θ = [ u 1 T , L , u d T , λ 1 , L , λ d , τ ] T . ]]>假设有PDF参数簇
f ( X | θ ( k ) ) = Π t = 1 n 1 π m | R ( k ) | exp ( - x t H [ R ( k ) ] - 1 x t ) - - - ( 6 ) ]]>
其中对应第k个候选模型。令H_k表示为假设信源个 数是k，
根据贝叶斯法则，可以得到
f ( H k | X ) = f ( X | H k ) f ( H k ) f ( x ) · - - - ( 7 ) ]]>
假设服从均匀分布，且f(X)独立于k，k值不影响公 式(7)的最小化，因此，由公式(7)可得
max k ∈ [ 0 , m ‾ - 1 ] f ( H k | X ) = max k ∈ [ 0 , m ‾ - 1 ] f ( X | H k ) - - - ( 8 ) ]]>
公式(8)说明H_k假设下最大化检测概率等效于最大化后验(MAP)概率。后验 概率(APP)计算为

其中f(X，θ^(k))表示X和θ^(k)的联合PDF，f(θ^(k))表示θ^(k)的先验PDF，是θ^(k)的最大似然估计，ν_k是θ^(k)的长度

是Hessian矩阵。通过计算的数学期望得到费雪信息矩阵
J = - E [ ∂ 2 log f ( X | θ ( k ) ) ∂ θ ( k ) ∂ ( θ ( k ) ) H ] · - - - ( 11 ) ]]>
由公式(9)得

近似符号是由于logf(θ^(k))和(ν_k/2)log2π独立于n，且m固定n→∞时 这里O(1)表示n→∞时该值趋于一个常数。接下来忽略独立于k 的项并令ν_k＝k(2m-k)，经典BIC方法可写成

针对k值最小化公式(13)即可以估计信源数目。值得注意的是，公式(13)的 准则也可由基于MDL的原理得到。
由于m，n→∞和m/n→c，观测信息矩阵不仅仅依赖于n而且依赖于m。 因此公式(12)近似不再有效，公式(13)中经典BIC方法的性能将大大下降， 尤其是在快拍与天线数量相当时。为了绕过这个问题，重新计算m，n→∞和 m/n→c情况下的LLF(log-likelihood function，对数形式的似然函数)和 PF，则可以得到能准确检测出大规模MIMO中信源数目的改良BIC方法。
为了正确计算广义渐近条件下，m，n→∞和m/n→c时，信源数估计的 APP，我们首先确定向量参数θ^(k)的ML估计，附录1证明了广义渐近条件 下θ^(k)的ML估计与经典渐近条件下估计相同

是m，n→∞和m/n→c时θ^(k)的ML估计，其中另一方面， 附录2表明随着m，n→∞和m/n→c，log-APP可以计算为

为了确定采样特征向量的近似行为，描述由特征向量确定的子空间是没有 意义的，因为其维数随m→∞增加到无穷。相反，确定特征投影矩阵对应 的二次函数是可行的。同样的，这里不单独讨论向量参数θ^(k)，其维数将随 m→∞增加到无穷，转而讨论公式(15)中函数的最大化，其最大化能有 效提高信源数目的检测概率。
已知代表向量参数θ^(k)的先验PDF，n→∞时有界。远 小于因为后者随m→∞或n→∞无上界。另一方面，随m→∞和 n→∞，ν_klogπ也远小于因而由公式(15)得到

当对k最小化时，如果m→∞和n→∞情况下能正确计算出和该式能准确估计出信源数目。因为是θ^(k)在m→∞和n→∞情况下的ML估计，所以LLF可以计算为

另一方面，由(B.4)可得的行列式

其中
Q = Q 11 Q 12 Q 21 Q 22 - - - ( 19 ) ]]>
Q₁₁，Q₁₁，Q₁₁和Q₁₁定义在(B.8)中。利用分块矩阵求行列式的公式可得
| Q | = | Q 11 | × | Q 22 - Q 21 Q 11 - 1 Q 12 | . - - - ( 20 ) ]]>
将(B.8)带入(20)得
| Q | = ( 2 n ) mk n k ( Π i = 1 k 1 l i ) m - k + 2 τ ^ k ( m - k ) k - - - ( 21 ) ]]>
然后将该式带入(18)，推导出

两边同时取对数，可以得到

由于m，n→∞，而假设信源数k是一个固定值。因此可得[log((m-k)n)]/m→0， (klogn)/m→0，(k-2)/m→0和(k²+2)/m→0，m，n→∞，m/n→c，(23)式 近似为

将(17)和(24)带入(16)，本发明所提出的改良BIC为
BIC ( k ) = 2 n ( m - k ) log a ( k ) g ( k ) + P ( k , m , n ) - - - ( 25 ) ]]>
其中分别是算术平均值和几何平均值。信 源数估计为
d ^ = arg min k = 0 , L , m ‾ - 1 BIC ( k ) . - - - ( 26 ) ]]>
已知中是最大假定信源数，其中m，n→∞，m/n→c。一对互 转置矩阵的特征值相同，因而天线数m和采样数n可以互换。这意味着m＞n 时，本发明的改良BIC方法应用于n个非零特征值时，可交换m和n。附 图1是本发明的在大规模自适应天线下基于贝叶斯信息准则的信源数估计 方法的流程图。
已知l₁，L，l_k是信号总体特征值的ML估计，是对应的信噪 比(SNR)。说明P(k，m，n)不仅与快拍数n有关，也与天线数m和SNR有关。 这就是说，PF项P(k，m，n)比在标准BIC中承载了更多信息，在m，n→∞和 m/n→c情况下计算结果更加精确。
为证明本发明方法的优势，比较了典型的ICT，BIC，基于线性收缩的 MDL(LS-MDL)，EEF，RMT-AIC和BN-AIC，其中BN-AIC的用户自定义 参数C设置为2。考虑一均匀线阵，阵列间距为半波长，分别接收窄带， 等功率平稳高斯信号。附图2-5中，″proposed″代表本发明所提出的方法。
实验1 正确检测概率与SNR的关系
在这个仿真中，考虑小样本n＝60，天线数是15。由附图2(a)看出，本 发明所提出的改良BIC方法的检测概率优于ITC。如附图2(b)所示，当快 拍数足够大，n＝150时，本发明的方法与现有ITC方法的差异变小。在大 样本情况下，本发明的方法性能与EEF相当，仍然超出LS-MDL和 RMT-AIC方法约0.5dB。而且，相比标准的BIC方法，本发明的方法的检 测性能明显提高了。
实验2 正确检测概率与DOA间距的关系。
附图3描述了小样本和大样本情况下正确检测概率与DOA间距的关 系。考虑三个不同信号的DOA，天线数为15。由附图3(a) 和附图3(b)可知，本发明的方法相比ITC方法，信源检测结果更精确， 尤其是在样本数较小时。很明显，这是由于标准BIC的代价函数项尚未优 化，低估信源数的概率占主导地位。本发明的方法给出了更小的代价函数， 可以减小低度拟合的概率，显著增强了检测性能。
实验3 正确检测概率与信源数的关系
该实验显示了大天线数与小快拍情况下的检测性能。考虑m＝40和 n＝80，附图4(a)和附图4(b)分别画出小信源数d＝3和大信源数d＝8情况下 的正确检测概率。由附图4(a)可以看出，本发明的方法明显优于ITC方法。 尽管BN-AIC方法在SNR足够大的情况下能够得到正确检测概率，但是结 果远远不如本发明的方法结果精确。附图4(b)表明，d＝8条件下，本发明 的方法优于ITC方法，在大信源数情况下，EEF方法不能获得正确检测概 率，除非SNR变得足够大。而且，由附图4(b)可知，LS-MDL方法比BIC 和RMT-AIC方法更有效，但次于本发明的方法。另外，本发明的方法的 检测性能显著超过了BN-AIC方法。
实验4 正确检测概率与快拍数的关系
该实验中m和n以相同速度增长且m/n＝0.3和m/n＝0.5两种情况下比较 算法性能。天线数和快拍数无线增长，但信源数不变，令m 和n以相同速度c＝m/n增长。附图5(a)表明EEF和LS-MDL算法比标准 BIC检测方法的结果更加精确，然而这些算法都不能在快拍数变得很大时 得到信源数的一致估计。比较ITC方法，本发明的方法能够得到更加精确 的结果。当m/n＝0.5，d＝8，本发明的方法的正确检测概率相比ITC算法 更快收敛到1，如附图5(b)所示。
A.附录1：θ^(k)在广义渐近条件下的ML估计
当信源数为k，令R^(k)＝UΛU^H，分别是R^(k)和的特征值分解 (EVDs)。其中，Λ＝diag(λ₁，L，λ_k，τ，L，τ)，U＝[u₁，L，u_m]，L＝diag(l₁，L，l_m)， E＝[e₁，L，e_m]，u_i和e_i，i＝1，L，m分别是总体特征值λ_i对应的总体特征向量和 采样特征值l_i对应的采样特征向量。因此，计算LLF

其中tr[·]代表求迹，G＝E^HU，因为G是正交的，所以得到下面不等式
tr ( Λ - 1 G H LG ) ≥ Σ i = 1 m l i λ i - - - ( A . 2 ) ]]>
(A.2)满足时，G＝I_m，即U＝E。由(A.1)和(A.2)可知，e_i，i＝1，L，m是u_i的 最大似然估计，即i＝1，L，k。将上述ML估计带入(A.1)，得到向量 参数简化的

对最大化λ₁，L，λ_k，τ的ML估计为

τ ^ k = 1 m - k Σ i = k + 1 m l i · - - - ( A . 5 ) ]]>
因此，得到广义渐近条件下ML估计该值与经典 渐近情况相同。
B.附录1：公式(15)的推导
公式(B.1)是的泰勒级数展开，其中是(10) 中定义的Hessian矩阵，为了简化删掉上标(·)^(k)。

下面证明随m，n→∞和m/n→c，(B.1)中零阶项远大于二阶项。由于 是θ在广义渐近条件下的ML估计。(B.1)中零阶项 可以表示为

为了确定(B.1)中的二阶项，我们要求-log f(X|θ)对θ的二阶偏导：

Hessian矩阵计算如下

如果且λ_i的重数为1，m，n→∞，m/n→c，则
l i = λ i + λ i τ λ i - τ c + O ( 1 n ) - - - ( B . 5 a ) ]]>
τ ^ k = τ + O ( 1 n ) · - - - ( B . 5 b ) ]]>
另一方面，在相同条件下，最大采样特征向量和总体特征向量的内积 几乎然收敛到一确定值

接下来，令Δθ@[ò^T，ν^T，ε]^T，其中
ò＝[(u₁-e₁)^T，L，(u_k-e_k)^T]^T    (B.6a)
ν = - τc [ λ 1 λ 1 - τ , L , λ k λ k - τ ] T + O ( 1 n ) - - - ( B . 6 b ) ]]>
ϵ = O ( 1 n ) - - - ( B . 6 c ) ]]>
(B.1)中二阶项表达式为

其中

Q 12 = blkdiag ( - 2 n e 1 l 1 , L , - 2 n e k l k ) - - - ( B . 8 b ) ]]>
Q 22 = diag ( n l 1 2 , L , n l k 2 ) - - - ( B . 8 c ) ]]>
β = n ( m - k ) τ ^ 2 - - - ( B . 8 d ) ]]>
这里，blkdiag(·)表示块对角矩阵。把(B.6)和(B.8)带入(B.7) 计算出特别地

其中和因为e₁，L，e_m和u₁，L，u_m张成了相同的观测空间，因此断言，对u_i(i＝1，L，m)，存在非零集合 {α_i1，L，α_im}，使得
u_i＝α_i1e₁+L+α_ime_m    (B.10)
这意味着
u i H u i = | α i 1 | 2 + L + | α im | 2 = 1 - - - ( B . 11 ) ]]>
其中|α_ij|表示α_ij的绝对值，容易得到



将(B.12)带入(B.9)得

(B.7)的二阶项和三阶项为

另外，容易得(B.7)的后两项
ν H Q 22 ν = n Σ i = 1 k ( λ i τc ) 2 l i 2 ( λ i - τ ) 2 - O ( n ) - - - ( B . 13 c ) ]]>
ϵ 2 β = ( m - k ) O ( 1 n ) = O ( 1 ) . - - - ( B . 13 d ) ]]>
最后，将(B.13)带入(B.7)，得到

利用l₁≥L≥l_m和可以推断出

然而，从(B.2)知，随m，n→∞和m/n→c，有界，忽略(B.1)中 高阶项可得

另一方面，假设θ的先验PDF围绕可得将该结果和(B.16) 带入(9a)

(B.17)取对数最终推导出(15)。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说 明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术 领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若 干简单推演或替换，都应当视为属于本发明的保护范围。