一种基于阵列信号处理的波达方向估计方法

1.本发明属于波达方向定位技术领域，尤其涉及一种基于阵列信号处理的波达方向估计方法。


背景技术：

2.阵列信号处理作为信号处理领域一个重要的研究方向，近年来得到了广泛的关注与迅速的发展。该项技术通过在三维空间中按一定流型安放的信号接收传感器阵列来获得信号的离散观测数据，并利用不同传感器接收数据的差异，设计相应的算法以获得信号的各项参数。近年来，阵列信号处理已被广泛应用于声呐通信、射电天文、生物医学、定位导航等众多民用或军用领域。
3.阵列信号处理的研究主要包括自适应空域滤波和空间谱估计两个方面。其中，空间谱估计也称为波达方向(direction of arrival,doa)估计。doa估计问题是根据阵列参考天线接收到的各个信号信息确定同时位于空间某片区域内多个信号的空间位置，也即是各个信号相对于阵列参考天线的方向角。doa估计的分辨率和天线阵列的长度有关，当天线阵列的长度确定之后，其对应的分辨率也就被确定，称其为瑞利限。对于瑞利限而言，存在能够突破瑞利限限制的方法，这类方法称为超瑞利限方法或超分辨方法。
4.早期的doa估计研究只研究空间来波信号的一维到达角度，接收阵列通常采用均匀线阵。最早提出的是常规波束形成法。随后，许多能够突破“瑞利限”且可以实现更高分辨率的估计方法先后被许多人提出来。其中以多重信号分类(music)算法最为出名。在1989 年,旋转不变技术(esprit)算法被roy r等人提出来。20世纪80年代，出现了诸如最大似然(ml)算法、子空间拟合(sf)算法等子空间类拟合算法。这些算法与music算法和 esprit等算法相比，估计的角度值更加准确，但是运算复杂度也更大。
5.近年来，压缩感知(cs)理论的发展也为阵列信号处理中doa估计提供了新思路。cs 理论跳出奈奎斯特定理的限制，当以低于奈奎斯特频率的采样频率对信号采样时，其能从维数不多的采样信号中完整恢复出原有信号。cs理论被应用在在空间谱估计中时，接收信号在整个空域来说可以看成是以低于奈奎斯特频率进行采样出的信号，然后通过求解最优解的方式即能找出信源到达角度。与传统的doa估计方法相比基于压缩感知理论的doa (cs-doa)估计具有众多优点，包括采样率低使得采样时间更短，抗干扰能力强等优势。目前，cs-doa估计仍处于不断探索的阶段，但其理论研究价值和实际工程意义不可估量。


技术实现要素：

6.本发明目的在于提供一种基于阵列信号处理的波达方向估计方法,通过天线阵列获得来波采样数据，由循环平稳算法计算循环相关矩阵，接着利用奇异值分解获得对应的信号子空间，构建稀疏重构模型，然后基于旋转不变子空间信息获得每次迭代残差对应的字典向量，将字典向量匹配后加入支撑集中计算更新残差，不断迭代直至计算出所有来波方向，以解决不能满足室内定位对精度和时延所提出要求的技术问题。
7.为解决上述技术问题，本发明的具体技术方案如下：
8.一种基于阵列信号处理的波达方向估计方法，包括以下步骤：
9.步骤1、通过天线阵列获得来波采样数据；
10.步骤2、由循环平稳算法计算循环互相关矩阵；
11.步骤3、利用奇异值分解获得对应的信号子空间，构建稀疏重构模型；
12.步骤4、基于旋转不变子空间信息获得每次迭代残差对应的字典向量；
13.步骤5、将字典向量匹配后加入支撑集中计算更新残差；
14.步骤6、不断迭代直至计算出所有来波方向，重复步骤4和步骤5，不断迭代直至计算出所有来波方向。
15.进一步的，步骤1具体包括以下步骤：采用均匀线阵对信号进行接收，对于相隔d 共有n个阵元的均匀线阵，f0为信号频率，λ0为信号波长；当有m个循环频率为α的信号源入射时，其中m≤n，获得第k个阵元在第t个快拍处接收的信号sk(t)，该信号来波方向θi满足-π/2≤θi≤π/2，其中i＝1,2，
…
，m，表示共有m个信号源；获得阵列采样的数据s(t)＝a(θ)x(t)+n(t)，t＝1,2,......,u，其中a(θ)表示接收信号角度矩阵，x(t)表示基础信号矩阵，n(t)表示噪声矩阵，t表示阵列的第t个快拍的数据，设置共接收了u 个快拍数据。
16.进一步的，步骤2的计算方法包括以下步骤：
17.步骤2.1、计算第k个阵元和第i个信号源的时延τ
ki
；
18.步骤2.2、由循环互相关的定义及阵元接收数据可得阵元s
p
和sq上接收数据的循环互相关函数
19.步骤2.3、由中p,q＝1,2,.,n，对n个不同阵元求解循环互相关函数，获得阵列接收信号的循环互相关矩阵其中p,q表示第p和第q个阵元
20.进一步的，所述步骤2.1中，
21.步骤2.2中，
22.其中表示第i个信号源基础信号的循环自相关，τ
pi
表示第p个阵元和第i个信号源的时延，τ
qi
表示第q个阵元和第i个信号源的时延，且对于等距均匀线阵而言,有下式成立：
[0023][0024]
由于噪声循环频率不为α，其在α频率上的相关系数为0，得到循环互相关函数为：
[0025][0026]
其中s
p
(t)表示第p个阵列采样的数据，sq(t)表示第q个阵列采样的数据；
[0027]
步骤2.3中阵列接收信号的循环互相关矩阵计算为：
[0028][0029]
进一步的，步骤3包括以下步骤：
[0030]
对步骤2所得循环互相关矩阵进行奇异值分解，获得酉矩阵u和v；将奇异值由大到小按顺序排列，根据奇异值大小将酉矩阵u和v分成信号子空间和噪声子空间两部分；从m个最大奇异值的特征向量中张成信号子空间us和vs；最终得到新的信号输出方程rs。
[0031]
进一步的，所述步骤3中，循环互相关矩阵的第i列像稀疏恢复模型一样重构：其中i＝1,2,...,n；和ri均表示的第i列，为字典矩阵，bi是字典下ri的稀疏系数向量；表示噪声值，ei是一个n
×
1列向量，第i个元素为1，其余为0；
[0032]
根据的所有列将其更改为矩阵形式为：其中 b＝[b1,b2,...,bn]并且bi(i＝1,2,..,n)具有相同的稀疏结构，即每个理想值的非零元素应该出现在b的同一行中，填充0的行数相同，in则是n
×
n的单位矩阵；
[0033]
其中是一个n
×
n的矩阵；的奇异值分解表示为：其中 u和v均为酉矩阵，σ由奇异值构成；其中因此有m个较大奇异值；将奇异值由大到小按顺序排列：λ1≥λ2≥...≥λm≥λ
m+1
≥...≥λn；后n-m个奇异值为噪声和干扰导致的；根据奇异值大小将酉矩阵u和v分成两部分：其中σs由前m大的奇异值构成，σn由剩余的奇异值构成；
[0034]
其中us和vs表示从m个最大奇异值的特征向量中张成的信号子空间，σ对角线中的其余n-m个奇异值构成噪声子空间un和vs；可以得到一个新的信号输出方程为：
[0035][0036]
其中为噪声干扰误差项；δr是由与的近似引起的误差；是信号子空间的组成部分。
[0037]
进一步的，步骤4计算方法包括如下步骤：
[0038]
步骤4.1、初始化状态支撑集s0＝φ，残差c0＝rs，稀疏表示矩阵
[0039]
步骤4.2、设此时为第k次迭代：支撑集为sk，残差为ck，稀疏表示矩阵为第 k次迭代的前p个最显著解为θi,i＝1,2,...,p；得a(θ)
·
t＝ck，其中t为旋转矩阵；设
[0040][0041]
其中a(θ)表示接收信号角度矩阵，a1(θ)表示a(θ)的前n-1行矩阵，a2(θ)表示a(θ)
的后n-1行矩阵，表示ck的前n-1行矩阵，分别表示ck的后n-1行矩阵；
[0042]
设过渡矩阵ψ，其中a2(θ)＝a1(θ)
·
ψ；计算得到由于t-1
ψt与ψ相似且ψ为对角矩阵，即可解得ψ；从而计算得到第k次迭代的字典向量αi,i＝1,2,...,p。
[0043]
进一步的，步骤5包括以下步骤：设第k次迭代的所选择的向量角标为ik，取将ik添加至sk，更新支撑集和稀疏表示矩阵；从而更新正交匹配追踪的残差迭代方程，获得新的残差解c
k+1
。
[0044]
进一步的，所述步骤5中，更新支撑集和稀疏表示矩阵：s
k+1
＝sk∪{ik}，
[0045]
正交匹配追踪的残差迭代方程表示为：
[0046][0047][0048]
其中下标k表示第k次迭代，b表示残差迭代方程的未知数，bk表示第k次迭代的最小二乘解，ck表示第k次迭代的残差信号矩阵，rs表示步骤3所获得的信号输出方程，表示第k次迭代时已经选定的k个原子所构成的新的字典矩阵，表示的伪逆矩阵，f指二范数；由上式可得对应残差c
k+1
，进行迭代。
[0049]
本发明的一种基于阵列信号处理的波达方向估计方法，具有以下优点：
[0050]
本发明基于循环平稳算法，旋转不变子空间信息以及正交匹配追踪算法，得到一种基于阵列信号处理的波达方向估计方法。利用循环平稳算法，将时间刻度上的循环平稳性与空间刻度上的阵列信号相结合，使得算法的抗噪性能增强。所述方法通过旋转不变子空间信息与基于压缩感知理论的正交匹配追踪算法的结合，突破正交匹配追踪算法存在的瑞利极限，有效提高多物体定位分辨力并降低计算复杂度，并进一步提高算法的抗干扰性。达到了室内定位场景高精度，低时延的要求，改善用户体验。
附图说明
[0051]
图1为本发明的基于阵列信号处理的波达方向估计方法的场景示意图；
[0052]
图2为本发明的基于阵列信号处理的波达方向估计方法流程图；
[0053]
图3为本发明的基于阵列信号处理的波达方向的cyclic-正交匹配追踪算法流程图；
[0054]
图4为本发明的基于阵列信号处理的波达方向的基于旋转不变子空间信息的精简字典正交匹配追踪算法流程图。
具体实施方式
[0055]
为了更好地了解本发明的目的、结构及功能，下面结合附图，对本发明一种基于阵列信号处理的波达方向估计方法做进一步详细的描述。
[0056]
如图1所示，本发明的研究基于室内定位场景。其中绑定在所需定位物体上的发信器用于发送信号，多个不同位置的天线阵列用于采样接收信号，处理器用于根据采样信号进行操作以执行所述的方法，将所得角度输入pc处理端后，根据多个角度数据获得高精度定位信息。
[0057]
如图2所示，本发明设计了一种基于阵列信号处理的蓝牙波达方向的方法，通过该方法对信号进行处理，获得对应的来波方向，用于实现高精度，低时延的室内定位。实际应用当中，具体包括如下步骤。
[0058]
步骤1、通过天线阵列获得来波采样数据。我们采用均匀线阵对信号进行接收，对于相隔d共有n个阵元的均匀线阵，当仅有一个信号源入射，忽略噪声时，其第一个阵元接收的信号为：其中a(t)表示接收信号的幅度，e为自然对数函数的底数， j为虚数，t为对应时间。相邻阵元间时延为其中c为光速，θ为对应来波方向。在窄带信号中，第二个阵元接收信号为：
[0059][0060]
其中，f0为信号频率，λ0为信号波长。
[0061]
当有m个循环频率为α的信号源入射时(m≤n)，第k个阵元在第t个快拍处接收的信号为：
[0062]
其中其中表示第k个阵元在所接受的第i个来波方向的信号，ai(t)表示接收信号的幅度，表示第 k个阵元在接收第i个来波方向时的时延，nk(t)为第k个阵元接收的噪声。其来波方向θi满足-π/2≤θi≤π/2，其中i＝1,2，
…
，m。令则阵列在第t个快拍处采样的数据可以表示为：
[0063][0064]
即：s(t)＝a(θ)x(t)+n(t)，t＝1,2,......,u。其中，u是快拍的数量， a(θ)＝[a(θ1),...,a(θm)]，a(θi)是第i个源的方向向量。
[0065]
在步骤2和步骤3中，阐述了本发明设计的基于阵列信号处理的蓝牙波达方向 cyclic-正交匹配追踪算法，如图3所示。
[0066]
步骤2、由循环平稳算法计算循环相关矩阵。本算法采用了cyclic算法进行计算。设时延τ
ki
，其中k表示第k个阵元，i表示第i个信号源，由循环互相关的定义及阵元接收数据可得阵元s
p
和sq上接收数据的循环互相关函数为
[0067][0068]
其中p＝1,2,
…
,n,q＝1,2,
…
,n,表示第i个信号源基础信号的循环自相关，τ
pi
表示第p个阵元和第i个信号源的时延，τ
qi
表示第q个阵元和第i个信号源的时延。且对于等距均匀线阵而言,有下式成立：
[0069][0070]
由于噪声循环频率不为α，其在α频率上的相关系数为0，得到循环互相关函数为：
[0071][0072]
其中s
p
(t)和sq(t)分别表示第p和q个阵列采样的数据。
[0073]
由于中p,q＝1,2,.,n，则阵列接收信号的循环相关矩阵可以计算为：
[0074][0075]
步骤3、利用奇异值分解获得对应的信号子空间，构建稀疏重构模型，循环相关矩阵的第i列像稀疏恢复模型一样重构：其中i＝1,2,...,n；和ri均表示的第i列，为字典矩阵，bi是字典下ri的稀疏系数向量；表示噪声值，ei是一个n
×
1列向量，第i个元素为1，其余为0；
[0076]
根据的所有列将其更改为矩阵形式为：其中 b＝[b1,b2,...,bn]并且bi(i＝1,2,..,n)具有相同的稀疏结构，即每个理想值的非零元素应该出现在b的同一行中，填充0的行数相同，in则是n
×
n的单位矩阵；
[0077]
其中是一个n
×
n的矩阵；的奇异值分解可以表示为：其中u和v均为酉矩阵，σ由奇异值构成；其中因此有m个较大奇异值；将奇异值由大到小按顺序排列：λ1≥λ2≥...≥λm≥λ
m+1
≥...≥λn；后n-m个奇异值为噪声和干扰导致的；根据奇异值大小将酉矩阵u和v分成两部分：其中σs、σn分别由前m大的奇异值和剩余的奇异值构成。
[0078]
其中us(vs)表示从m个最大奇异值的特征向量中张成的信号子空间，σ对角线中的其余n-m个奇异值构成噪声子空间un(vs)；结合上方两式中所示的关系，可以得到一个新的信号输出方程为：
[0079]
[0080]
其中为噪声干扰误差项；δr是由与的近似引起的误差；是信号子空间的组成部分。
[0081]
在步骤4，步骤5和步骤6中，阐述了本发明设计的基于阵列信号处理的蓝牙波达方向的基于旋转不变子空间信息的精简字典正交匹配追踪算法，如图4所示。
[0082]
步骤4、基于旋转不变子空间信息获得每次迭代残差对应的字典向量。
[0083]
设初始化状态支撑集s0＝φ，残差c0＝rs，稀疏表示矩阵
[0084]
设此时为第k次迭代：支撑集为sk，残差为ck，稀疏表示矩阵为第k次迭代的前p个最显著解为θi,i＝1,2,...,p；
[0085]
令span(a(θ))、span(us)分别表示由a(θ)和us所张成的空间，由于 span(a(θ))＝span(us)，由步骤3中信号输出方程rs可得 span(a(θ))＝span(rs)＝span(c0)；
[0086]
由残差迭代方程可知span(ck)＝span(rs)；
[0087]
由上述两式可得span(ck)＝span(a(θ))，所以ck与a(θ)在同一张成的子空间中；
[0088]
由上式可得a(θ)
·
t＝ck，其中t为旋转矩阵。设
[0089][0090]
其中a1(θ)和a2(θ)均有n-1行；考虑信号源数p，a1(θ),a2(θ)∈c
n-1
×
p
，得到 a2(θ)＝a1(θ)
·
ψ，其中
[0091]
同理，设
[0092][0093]
由上述三式可得可得
[0094]
由上式可得
[0095]
由上式可得
[0096]
由于t-1
ψt与ψ相似，即可解得ψ；
[0097]
可得字典向量：
[0098]
其中i＝1,2,...,p。
[0099]
步骤5、将字典向量匹配后加入支撑集中计算更新残差。设第k次迭代的所选择的向量角标为ik，取将ik添加至sk，更新支撑集和稀疏表示矩阵：
[0100]sk+1
＝sk∪{ik}，
[0101]
omp(正交匹配追踪)的残差迭代方程表示为：
[0102][0103][0104]
其中下标k表示第k次迭代，b表示残差迭代方程的未知数，bk表示第k次迭代的最小二乘解，ck表示第k次迭代的残差信号矩阵，rs表示步骤3所获得的信号输出方程，表示第k次迭代时已经选定的k个原子所构成的新的字典矩阵，表示的伪逆矩阵。由上式可得对应残差c
k+1
，进行迭代。
[0105]
步骤6、不断迭代直至计算出所有来波方向。重复步骤4和步骤5，不断迭代直至计算出所有来波方向。
[0106]
可以理解，本发明是通过一些实施例进行描述的，本领域技术人员知悉的，在不脱离本发明的精神和范围的情况下，可以对这些特征和实施例进行各种改变或等效替换。另外，在本发明的教导下，可以对这些特征和实施例进行修改以适应具体的情况及材料而不会脱离本发明的精神和范围。因此，本发明不受此处所公开的具体实施例的限制，所有落入本技术的权利要求范围内的实施例都属于本发明所保护的范围内。