一种存在阵列互藕的一比特采样DOA估计方法
processing(icassp).ieee,2017》中，该方法被扩展到稀疏阵，利用更多的虚拟阵元提高了doa估计性能。此外，稀疏贝叶斯学习，联合稀疏表示等方法在一比特doa估计中也有应用。
5.虽然互耦或一比特采样的问题已经在很多文献中得到了广泛的讨论，但在这两个问题同时存在的特定情况下的doa估计问题还没有得到解决。因此，研究如何在一比特采样下进行稳健的doa估计是十分必要的。
技术实现要素：
6.本发明的目的在于，提出一种存在阵元间互耦效应时的基于一比特采样的稳健doa估计方法，用于避免高精度模数转换器带来的高硬件成本和高系统功耗。该方法没有损失阵列孔径，且不需要阵列校准，实际操作起来比较简单。
7.本发明提供了一种存在阵列互藕的一比特采样doa估计方法，所述方法包括：
8.步骤1，设置天线阵列，所述天线阵列为采用m个阵元组成的均匀线性阵列，其中阵元间距为d＝λ/2，λ表示信号波长，定义天线阵列的接收信号数k，阵列的未量化观测数据矢量为x(t)，阵列的一比特观测数据矢量为y(t)；
9.步骤2：计算归一化协方差矩阵：
10.(2a)计算观测数据y(t)在快拍数为n的协方差矩阵ry；
11.(2b)计算归一化协方差矩阵
[0012][0013]
其中，矩阵是一个对角矩阵，中间m-2p+2个非零元素的取值相同，定义为γ，p表示允许忽略互藕效应的阵元间隔数，即相距p个阵元间隔的两个阵元之间的互耦效应可以忽略；和分别代表取复数的实部和虚部，j表示虚数单位，c表示互藕矩阵，且a表示k个信号的阵列流形矩阵，表示噪声功率，rs表示信源的协方差矩阵。
[0014]
步骤3，矩阵变换处理：
[0015]
(3a)将互耦矩阵和导向矢量的乘积转化为ca(θ)＝t(θ)α；
[0016]
其中t(θ)＝blkdiag{t1,t2,t3}，t1＝diag{1,β(θ),...,β
p-2
(θ)}，t2＝[β
p-1
(θ),...,β
m-p
(θ)]
t
；
[0017]
t3＝diag{β
m-p+1
(θ),...,β
m-1
(θ)}，β(θ)＝e
j2πdsinθλ
，α＝[μ1,...,μ
p-1
,τ(θ),α1,...,α
p-1
]
t
；；
[0018]
(3b)根据归一化矩阵g-1/2
和矩阵t(θ)的结构和元素大小，将二者乘积转化为g-1/2
t(θ)＝t(θ)d，得到g-1/2
ca(θ)＝g-1/2
t(θ)α＝t(θ)dα，其中d＝diag[ρ1,ρ2,...,ρ
p-1
,γ,ν1,ν2,...,ν
p-1
]是维度为2p-1的对角矩阵；
[0019]
步骤4，构建稀疏矢量：
[0020]
(4a)将归一化协方差矩阵列向量化，得到等效接收信号矢量r：
[0021][0022]
其中，其中，表示第k个信号的功率，vec(
·
)表示列向量化，
⊙
表示hadamard乘积。
[0023]
(4b)将等效接收信号矢量r变换表示为：
[0024][0025]
其中，η＝dα,其中，η＝dα,表示克罗内克积；
[0026]
(4c)将空域角度均匀划分为n个网格，记为其中n＞＞k；
[0027]
(4d)将等效接收信号矢量r重新表示为稀疏矢量：
[0028][0029]
其中，为字典矩阵，表示稀疏度为k(2p-1)2的块稀疏列向量；
[0030]
步骤5，获取准确噪声子空间：
[0031]
(5a)设置迭代次数i的初始值(例如i＝0)，设置收敛阈值ε(例如ε＝0.001)，初始化矩阵qi，例如qi＝im。im表示维度为m
×
m单位矩阵；
[0032]
(5b)将特征值分解，得到信号子空间e
is
和噪声子空间e
in
以及相应的特征值矩阵λ
is
和λ
in
；
[0033]
(5c)更新q
i+1
为其中表示以括号中矩阵的对角线元素构成的对角矩阵，并更新i＝i+1。
[0034]
(5d)重复步骤(5b)和(5c)直到满足
[0035]
步骤6，构造加权矩阵：
[0036]
(6a)根据空间网格角度分别计算其中e
in
为步骤5中最后一次迭代获得的噪声子空间，det(
·
)表示求行列式；
[0037]
(6b)构造加权矩阵其中w＝max{w1,...,wn}。
[0038]
步骤7，获取估计的doa：
[0039]
(7a)基于块稀疏矢量构建稀疏矢量其中，ho的第n个元素
表示的第l个元素；
[0040]
(7b)利用等效接收信号矢量构建优化问题：
[0041][0042]
其中，ξ为正则化参数，||
·
||1和||
·
||2分别表示l1和l2范数；
[0043]
对所构建的优化问题进行求解，求解出的稀疏矢量ho中最大的k个元素对应的角度即为估计的doa。
[0044]
本发明提供的技术方案至少带来如下有益效果：
[0045]
在本发明中，使用一比特采样，降低了系统成本和功耗以及数据传输存储量；无需预先的阵列校准即可实现稳健的doa估计，实现起来较为简单；没有损失阵列孔径，角度分辨率较高。
附图说明
[0046]
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。
[0047]
图1为本发明针对的均匀线性阵列结构设置示意图。
[0048]
图2为本发明实施方式仿真实验doa的均方根误差随snr变化关系示意图。
[0049]
图3为本发明实施方式仿真实验doa的均方根误差随快拍数变化关系示意图。
[0050]
图4为本发明实施方式仿真实验doa分辨概率随角度间隔变化关系示意图。
具体实施方式
[0051]
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。
[0052]
本发明针对的是阵列为均匀线性阵列，阵元数为m，阵元间距为d＝λ/2，其中λ为信号波长，具体阵列结构如图1所示。具体算法为：首先，利用阵列的一比特观测矢量计算出一比特协方差矩阵，再根据arcsine law恢复出归一化协方差矩阵。然后，为了对抗均匀线性阵列中未知的互耦效应，将互耦矩阵和导向矢量的乘积重新参数化，又利用了一个矩阵变换避免了归一化矩阵的影响。随后，通过网格划分构建了一个稀疏矢量。之后，利用迭代步骤解决了等效的非均匀噪声矩阵造成无法获得准确的噪声子空间的问题。接着，构造了一个加权矩阵用于增强稀疏恢复问题中稀疏矢量的稀疏性。最后，构建了一个单测量矢量的稀疏恢复问题进行doa估计。
[0053]
其具体实现步骤包括：
[0054]
步骤1：设置天线阵列：
[0055]
设置一个如图1所示的均匀线性阵列，共有m＝11个阵元，阵元间距d＝λ/2。假设共k个窄带平稳独立高斯信号入射到阵列上，信号的到达方向为{θ1,θ2,
…
,θk}，则阵列的未量化观测矢量在t时刻可以表示为：
[0056][0057]
其中，s(t)＝[s1(t),s2(t),...,sk(t)]
t
表示源信号波形矢量。表示每个元素均服从高斯分布的噪声矢量。代表k个信号的阵列流形矩阵，θk表示第k个信号的到达方向，其中代表第k个信号的导向矢量。
[0058]
若考虑阵元间互耦效应，将会使上述信号模型产生变化，特别地对于均匀线性阵列，互耦矩阵可以被描述为复对称带状托普利兹矩阵。而阵元间的距离越大互耦效应的影响就越小，当两个阵元之间的距离足够远时，可以近似为零。若假设当任意两个阵元被个阵元间隔分隔分开时，这两个阵元间的互耦效应是可以忽略的。那么，可以将互耦矩阵记为
[0059][0060]
这里且c1＝1＞|c2|＞
…
＞|c
p
|,c2～c
p
表示互耦矩阵c的第2至第p个元素。
[0061]
这种情况下，信号模型需要被改写为：
[0062][0063]
为了降低硬件成本和系统功耗，采用一比特adc对接收数据进行量化，则相应的阵列输出变为：
[0064]
y(t)＝signe(x(t))＝signe(cas(t)+n(t))
ꢀꢀꢀ
(9)
[0065]
其中，表示复值量化函数，这里sign(
·
)表示符号函数，和分别代表取复数的实部和虚部。
[0066]
步骤2：计算归一化协方差矩阵：
[0067]
首先x(t)的协方差矩阵可以被写为
[0068][0069]
其中，e[]表示数学期望，rs表示信源的协方差矩阵，rs＝e{s(t)sh(t)}，表示噪声功率，i表示单位矩阵。
[0070]
然后定义归一化协方差矩阵为：
[0071][0072]
其中g是一个对角矩阵，满足《g》
m,m
＝《r
x
》
m,m
。为了便于说明，定义空间相位因子根据等式(10)，g中与第m个天线(p-1＜m＜m-p+1)对应的对角元素可以表示为
[0073][0074]
由公式(12)可知g的第p-1个到第m-p+1个对角元素的值是相等的，这里记为γ，然后g就可以被写为其中ρ1,...,ρ
p-1
和ν1,...,ν
p-1
分别代表g的前p-1和后p-1个对角元素。
[0075]
观测数据y(t)的协方差矩阵ry的计算公式为：
[0076][0077]
其中，n为快拍数。y(tn)表示第tn个采样点的观测数据。
[0078]
当信源服从独立的高斯分布时，归一化协方差矩阵和一比特协方差矩阵ry之间的关系可以由arcsine law建立为
[0079][0080]
其中
[0081]
因此，可以用公式(14)重建归一化的协方差矩阵，表示为
[0082][0083]
其中
[0084]
步骤3：矩阵变换：
[0085]
为了对抗均匀线性阵列中未知的互耦效应，利用一个经典的变换将互耦矩阵和导向矢量的乘积转化为
[0086]
ca(θ)＝t(θ)α
ꢀꢀꢀ
(16)
[0087]
其中
[0088][0089]
是一个维度为m
×
(2p-1)的矩阵，而
[0090]
α＝[μ1,...,μ
p-1
,τ(θ),α1,...,α
p-1
]
t
ꢀꢀꢀ
(18)
[0091]
是一个长度为2p-1的向量，其中1的向量，其中
[0092]
基于以上变换，有g-1/2
ca(θ)＝g-1/2
t(θ)α。值得说明的是，矩阵左乘一个对角矩阵相当于矩阵的每一行分别与对角矩阵中相同行的对角元素相乘，矩阵右乘一个对角矩阵相当于矩阵的每一列分别与对角矩阵中相同列的对角元素相乘。而对角矩阵g-1/2
中有相同值的元素的行与t(θ)中间m-p+2行是对应的，又由于t(θ)的特殊结构，可以将t(θ)右乘一个对角矩阵转化为t(θ)左乘另一个对角矩阵。因此，有
[0093]
g-1/2
ca(θ)＝g-1/2
t(θ)α＝t(θ)dα
ꢀꢀꢀ
(19)
[0094]
其中
[0095]
步骤4：构建稀疏矢量：
[0096]
首先将归一化协方差矩阵列向量化，得到等效接收信号矢量；
[0097][0098]
其中为信源功率矢量，
[0099]
然后利用公式(19)中的矩阵变换将上式重新表示为
[0100][0101]
其中η＝dα,
[0102]
所以矩阵形式r可以被整理为：
[0103]
r＝bh+inꢀꢀꢀ
(22)
[0104]
其中
[0105]
随后将空间角度进行离散化，一共划分n个网格为其中n＞＞k,由此r可以被表示为稀疏矢量的形式：
[0106][0107]
其中是一个过完备字典矩阵，表示稀疏度为k(2p-1)2的列向量，每个块的维度为(2p-1)2。
[0108]
步骤5：获得准确噪声子空间：
[0109]
为了增强稀疏矢量的稀疏性，将构造一个加权矩阵，但公式(11)中的可以视为非均匀噪声矩阵，这将会影响特征值分解的结果，导致无法获得准确的子空间。为此，通过特征值分解，用k个大特征值及k个大特征向量构成的无噪声协方差矩阵对非均匀噪声矩阵进行更新，直至满足迭代终止条件，具体的迭代步骤如下：
[0110]
(5a)首先初始化迭代变量，设置i＝0，q0＝im，ε＝0.001。
[0111]
(5b)然后将特征值分解，有得到信号子空间和噪声子空间以及相应的特征值矩阵和(其中e
ks
和e
kn
分别为k个较大特征值和m-k个较小特征值对应的特征向量组成的矩阵，λ
ks
和λ
kn
为相应的特征值构成的对角矩阵)。
[0112]
(5c)更新q
i+1
为其中表示以括号中矩阵的对角线元素构成的对角矩阵，并设置i＝i+1。
[0113]
(5d)重复步骤(5b)和(5c)直到满足最后一次迭代中获得的e
in
即为准确的噪声子空间。
[0114]
步骤6：构造加权矩阵：
[0115]
考虑归一化协方差矩阵的无噪声部分
[0116][0117]
步骤5中最后一次迭代得到的e
in
即为的噪声子空间，然后根据子空间原理可知，对于真实的波达方向θk，矩阵是奇异的。
[0118]
然后根据空间网格角度分别计算分别计算det(
·
)表示求行列式。
[0119]
最后构造加权矩阵为其中w＝max{w1,...,wn}。
[0120]
步骤7：提出优化问题：
[0121]
知道信源的角度是由块稀疏矢量中非零元素的位置确定的，而中相同块的非零元素对应的角度也是相同的，根据这一特点构建稀疏矢量其中的第n个元素等于的第(n-1)(2p-1)2+1到第n(2p-1)2个元素的l2范数，即为范数，即为表示的第l个元素。
[0122]
利用加权矩阵，给出优化问题
[0123][0124]
其中，大的权重使得ho中对应的元素更有可能是0，小的权重使得ho中对应的元素得以保留，以此增强ho的稀疏性。
[0125]
对上述优化问题应用lasso模型，可以得到
[0126][0127]
其中ξ为正则化参数，||
·
||1和||
·
||2分别表示l1和l2范数，求解出的稀疏矢量ho中最大的k个元素对应的角度即为所估计的doa。
[0128]
为验证算法的性能，本实施例设置了三个仿真实验，仿真实验中将本发明所提方法与有互耦lasso,无互耦lasso，有互耦music,无互耦music以及中间子阵算法进行了比较。第一个仿真实验为几种算法在快拍数为2000的情况下，doa估计均方根误差随信噪比的变化关系。而第二个仿真实验为几种算法在信噪比为30db时，doa估计均方根误差随快拍数的变化关系。第三个仿真实验给出了信噪比为30db、快拍数为2000的条件下，本发明与中间子阵算法的对相邻信源的分辨概率随角度间隔的变化关系，实验的结果分别如图2、3、4所
示，图中的所提算法即指本发明实施例提供的存在阵列互藕的一比特采样doa估计方法。
[0129]
三个仿真实验的随机实验次数均为200，所用均匀线阵相同，阵元数同为11，阵列间距为d＝λ/2，互耦参数设置为p＝3，c2＝0.6864+0.4776i，c3＝0.4864+0.3376i，正则化参数为ξ＝0.5。前两个实验的入射的信号源数k＝2，到达角度分别为θ1＝10
°
和θ2＝55
°
。第三个实验同样考虑两个信号，入射角度分别为θ1＝10
°‑
δθ/2和θ2＝10
°
+δθ/2，其中δθ为角度间隔。本发明能在存在阵元间互耦效应时使用一比特采样数据进行doa估计，大大降低了系统的能量消耗和数据传输和存储的负担，并且利用了全部的阵列孔径，有着较高的doa估计性能和角度分辨率。
[0130]
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。
[0131]
以上所述的仅是本发明的一些实施方式。对于本领域的普通技术人员来说，在不脱离本发明创造构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。
