一种适用于分布式稀疏阵列的DOA估计方法

一种适用于分布式稀疏阵列的doa估计方法
技术领域
1.本发明属于doa估计算法技术领域。


背景技术：

2.稀疏阵列，如最小冗余阵列、互质阵列、嵌套阵列等，能够显著增加阵列的自由度来提高最大可分辨源信号的数量。分布式阵列具有很大的基线长度，有效增大了阵列的孔径，能够显著提高参数的估计精度，但其不能增加阵列的自由度，最大估计的源信号数目并没有增加。因此，分布式稀疏阵列结合了分布式阵列和稀疏阵列的优势，既能显著提高doa估计精度又能增加最大可分辨源信号的数量。在稀疏阵列的doa估计中，空间平滑广泛用于多重信号分类(music)算法和旋转不变子空间(esprit)算法中。然而，空间平滑算法复杂且计算量大，而且在稀疏阵列的doa估计中也不是必需的。因此，代替空间平滑，可以从差分协同阵列的接收数据中生成toeplitz矩阵(托普利兹矩阵)，估计效果与空间平滑算法相同。对于非均匀稀疏阵列，如嵌套阵列、互质阵列等，很多文章都详细介绍了toeplitz矩阵重构算法。与空间平滑算法相比，该方法大大降低了算法复杂度和计算复杂度。然而，对于分布式稀疏阵列，由于其包含多个子阵列，且子阵列之间的基线长度远大于子阵列的孔径，使用传统的toeplitz矩阵重构算法进行doa估计会产生模糊，无法估计出源信号的准确方位。而且，当分布式稀疏阵列的子阵数量和阵元天线数量相对较大时，采用空间平滑算法获得高精度的doa(波达方向)估计，由于其计算量巨大，既需要很长的运行时间，同时对计算机的硬件设备性能提出了更高的要求，这在很多实际工程中无法应用。


技术实现要素：

3.本发明目的是为了解决现有分布式稀疏阵列的doa估计算法存在计算量大，运算时间长，对硬件性能要求高的问题，提出了一种适用于分布式稀疏阵列的doa估计方法。
4.本发明所述一种适用于分布式稀疏阵列的doa估计方法，包括：
5.步骤一、对分布式稀疏阵列的阵元天线接收的信号进行采样，获得接收信号x(t)；
6.步骤二、计算接收信号x(t)的协方差矩阵r；
7.步骤三、对协方差矩阵r进行矢量化、去冗余和重新排列处理，获得协同阵列的接收信号，选择协同阵列连续部分的接收信号作为等效接收信号z1；
8.步骤四、根据所述等效接收信号z1的维度构造选择矩阵jq；
9.步骤五、利用选择矩阵jq选取等效接收信号z1中的数据构造toeplitz矩阵
10.步骤六、对toeplitz矩阵进行奇异值分解，得到信号子空间对应的特征向量，再采用多尺度esprit算法对所述信号子空间对应的特征向量进行求解，获取doa估计值。
11.进一步地，本发明中，步骤一中，接收信号x(t)表示为：
12.根据阵元天线间的互耦影响，以第一个稀疏子阵作为参考子阵，将第一个稀疏子阵的第一个天线作为参考天线，t时刻接收到的信号为：
[0013][0014]
其中，x(t)＝[x1(t),x2(t),
…
,xm·n(t)]
t
表示接收信号，xm(t)为第m个稀疏子阵的接收信号，m＝1,2,...,m，xi(t)表示第i个天线的接收信号，i＝1
,
2,...,m.n，上角标t表示矩阵转置，s(t)＝[s1(t),s2(t),
…
,sk(t)]
t
表示源信号，n(t)＝[n1(t),n2(t),
…
,nm·n(t)]
t
表示高斯白噪声，ni(t)表示第i个天线的加性高斯白噪声，是分布式稀疏阵列的阵列流形矩阵，as＝[a(θ1),a(θ2),
…
,a(θk)]表示参考稀疏子阵的阵列流形矩阵，θ1,θ2,
…
,θk分别表示第1至第k个源信号的波达方向；a(θ1)、a(θ2)分别为波达方向为θ1和θ2的稀疏子阵列的信号导向矢量，a(θk)表示波达方向为θk的稀疏子阵列的信号导向矢量，的稀疏子阵列的信号导向矢量，表示第n个天线的幅相值，n＝1,2,...,n，c为互耦矩阵，dn表示第n个天线与参考天线之间的间隔，λ表示载频波长，第m个稀疏子阵与参考稀疏子阵之间的相位差分矩阵表示为diag{
·
}表示对角矩阵，m表示稀疏子阵的总个数，n表示每个稀疏子阵列包含的天线个数，k表示待估计源信号的个数，d表示基线长度。
[0015]
进一步地，本发明中，步骤二中，接收信号x(t)的协方差矩阵r为：
[0016][0017]
其中，e{
·
}表示数学期望，上角标h表示共轭转置，表示噪声功率，rs为源信号的协方差矩阵，i表示单位矩阵。
[0018]
进一步地，本发明中，步骤三中对协方差矩阵r进行矢量化的方法为：
[0019]
应用khatri-rao积性质，将协方差矩阵r矢量化获得：
[0020][0021]
其中，vec(
·
)表示矢量化，上角标*表示取共轭，
⊙
表示khatri-rao(卡特里-拉奥)积，表示源信号的功率，分别表示第1个，第2个，
…
，第k个源信号的功率，z表示矢量化后得到的数据。
[0022]
进一步地，本发明中，步骤三中，等效接收信号z1：
[0023]
[0024]
其中，a1＝c0[a1(θ1),a1(θ2),
…
,a1(θk)]是分布式稀疏阵列协同阵列的中间子阵的阵列流形矩阵，c0为分布式稀疏阵列协同阵列的中间子阵的互耦矩阵，a1(θk)为波达方向为θk的均匀子阵列的信号导向矢量；l是协同阵列的子阵列的天线个数，表示中间元素为1其余为0的列向量，λ表示载频波长，z1可等效为包含2m-1个相同的均匀子阵列的分布式稀疏阵列的接收信号，该阵列称为分布式稀疏阵列的协同阵列，d表示天线之间的单位间距，等于半波长。
[0025]
进一步地，本发明中，步骤四中，构造选择矩阵jq为：
[0026]jq
＝[0
(2m-1)q
×
(q-1-q) blkdiag{y} 0
(2m-1)q
×q] (5)
[0027]
其中，0
(2m-1)q
×
(q-1-q)
表示(2m-1)q
×
(q-1-q)维的全零矩阵，q＝0,1,2
…
,q-1，y＝[i
(q)
,b1,b2,
…
,b
2m-2
]，b1＝b2＝
…
＝b
2m-2
＝[0q×
(q-1) i
(q)
]，jq∈{0,1}
(2m-1)q
×
(2m-1)l
，blkdiag{
·
}表示生成以矩阵块为对角线的对角矩阵，q为的秩，等于分布式稀疏阵列最大可检测源信号的数量。
[0028]
进一步地，本发明中，步骤五中，构造toeplitz矩阵为：
[0029][0030]
其中，jq为选择矩阵，q＝0,1,2
…
,q-1。
[0031]
对比本发明所述方法中采用的toeplitz矩阵重构算法与空间平滑算法的计算复杂度，如图1所示，包含m个子阵的分布式稀疏阵列，使用空间平滑算法为构造空间平滑矩阵首先计算q个协方差矩阵，需要(2m-1)2·
q3次复数乘法运算，然后执行约(2m-1)2·
q3次复数加法运算，最后运行(2m-1)2·
q2次除法。因此，空间平滑算法共需要执行o((2m-1)2·
q3)次复数乘法和o((2m-1)2·
q3)次复数加法运算。相比之下，本发明提出的适用于分布式稀疏阵列的新型toeplitz矩阵重构算法，可以直接使用选择矩阵jq选择接收信号z1中的数据构造toeplitz矩阵不需要任何算数运算，因此，本发明所述估计算法采用toeplitz矩阵重构算法，可以降低计算复杂度，特别是当子阵列数目和天线数目特别多时，可以大大减少运算量，提高计算速度。本发明所述算法既能获得与空间平滑算法相同的估计精度，又能降低算法的计算量，节省运行时间和计算成本。
附图说明
[0032]
图1是包含m个子阵的分布式稀疏阵列结构示意图；
[0033]
图2是包含m个子阵的分布式稀疏阵列的协同阵列示意图；
[0034]
图3是包含两个子阵的分布式超级嵌套阵列的协同阵列示意图；
[0035]
图4是本发明采用的toeplitz矩阵重构算法与现有空间平滑算法doa估计性能对比图；
[0036]
图5是本发明所述方法流程图。
具体实施方式
[0037]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完
整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0038]
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
[0039]
具体实施方式一：下面结合图1说明本实施方式，本实施方式所述一种适用于分布式稀疏阵列的doa估计方法，包括：
[0040]
步骤一、对分布式稀疏阵列的阵元天线接收的信号进行采样，获得接收信号x(t)；
[0041]
步骤二、计算接收信号x(t)的协方差矩阵r；
[0042]
步骤三、对协方差矩阵r进行矢量化、去冗余和重新排列处理，获得协同阵列的接收信号，选择协同阵列连续部分的接收信号作为等效接收信号z1；
[0043]
步骤四、根据所述等效接收信号z1的维度构造选择矩阵jq；
[0044]
步骤五、利用选择矩阵jq选取等效接收信号z1中的数据构造toeplitz矩阵
[0045]
步骤六、对toeplitz矩阵进行奇异值分解，得到信号子空间对应的特征向量，再采用多尺度esprit算法对所述信号子空间对应的特征向量进行求解，获取doa估计值。
[0046]
本实施方式中,协方差矩阵矢量化去冗余重新排列后的数据等效于如图2所示的协同阵列接收的单快拍数据。本发明提出了一种适用于分布式稀疏阵列的新型toeplitz矩阵重构算法，该算法能够获得与空间平滑算法相同的doa估计精度，但其计算复杂度可以大幅降低，同时降低了对硬件性能的要求，具有很好的实际工程应用前景。分布式稀疏阵列结构如图1所示，它包含m个相同的稀疏子阵列，用s1，s2，
…
，sm表示，所有子阵列均位于同一条直线上。每个稀疏子阵列包含n个阵元天线，总的天线数目为m
×
n。稀疏子阵列可以是嵌套阵列、超级嵌套阵列、互质阵列以及它们的扩展改进阵列。两个相邻稀疏子阵列之间的基线长度为d，其远大于稀疏子阵列的孔径。分布式稀疏阵列的天线位置集合表示为s＝s1∪s2∪...∪sm，∪表示集合的并集。
[0047]
进一步地，本实施方式中，结合图1进行说明，步骤一中，接收信号x(t)表示为：
[0048]
根据阵元天线间的互耦影响，以第一个稀疏子阵作为参考子阵，将第一个稀疏子阵的第一个天线作为参考天线，t时刻接收到的信号为：
[0049][0050]
其中，x(t)＝[x1(t),x2(t),
…
,xm·n(t)]
t
表示接收信号，xm(t)为第m个稀疏子阵的接收信号，m＝1,2,...,m，xi(t)表示第i个天线的接收信号，i＝1
,
2,...,m
.
n，上角标t表示矩阵转置，s(t)＝[s1(t),s2(t),
…
,sk(t)]
t
表示源信号，n(t)＝[n1(t),n2(t),
…
,nm·n(t)]
t
表示高斯白噪声，ni(t)表示第i个天线的加性高斯白噪声，是分布式稀疏阵列的阵列流形矩阵，as＝[a(θ1),a(θ2),
…
,a(θk)]表示参考稀疏子阵的阵列流形矩阵，θ1,
θ2,
…
,θk分别表示第1至第k个源信号的波达方向；a(θ1)、a(θ2)分别为波达方向为θ1和θ2的稀疏子阵列的信号导向矢量，a(θk)表示波达方向为θk的稀疏子阵列的信号导向矢量，的稀疏子阵列的信号导向矢量，表示第n个天线的幅相值，n＝1,2,...,n，c为互耦矩阵，dn表示第n个天线与参考天线之间的间隔，λ表示载频波长，第m个稀疏子阵与参考稀疏子阵之间的相位差分矩阵表示为diag{
·
}表示对角矩阵，m表示稀疏子阵的总个数，n表示每个稀疏子阵列包含的天线个数，k表示待估计源信号的个数，d表示基线长度。
[0051]
进一步地，本实施方式中，步骤二中，接收信号x(t)的协方差矩阵r为：
[0052][0053]
其中，e{
·
}表示数学期望，上角标h表示共轭转置，σ
n2
表示噪声功率，rs为源信号的协方差矩阵，i表示单位矩阵。
[0054]
进一步地，本实施方式中，步骤三中对协方差矩阵r进行矢量化的方法为：
[0055]
应用khatri-rao积性质，将协方差矩阵r矢量化获得：
[0056][0057]
其中，vec(
·
)表示矢量化，上角标*表示取共轭，
⊙
表示khatri-rao(卡特里-拉奥)积，表示源信号的功率，分别表示第1个，第2个，
…
，第k个源信号的功率，z表示矢量化后的得到的数据。
[0058]
进一步地，本实施方式中，步骤三中，等效接收信号z1：
[0059][0060]
其中，a1＝c0[a1(θ1),a1(θ2),
…
,a1(θk)]是分布式稀疏阵列协同阵列的中间子阵的阵列流形矩阵，如图2中的中间子阵的阵列流形矩阵，图中，ula1表示第1个均匀线性阵列，ula2表示第2个，...，ula
2m-1
表示第2m-1个均匀线性阵列，c0为分布式稀疏阵列协同阵列的中间子阵的互耦矩阵，a1(θk)为波达方向为θk的均匀子阵列的信号导向矢量；l是协同阵列的子阵列的天线个数，表示中间元素为1其余为0的列向量，λ表示载频波长，z1可等效为包含2m-1个相同的均匀子阵列的分布式稀疏阵列的接收信号，该阵列称为分布式稀疏阵列的协同阵列，d表示天线之间的单位间距，等于半波长。
[0061]
进一步地，本实施方式中，步骤四中，构造选择矩阵jq为：
[0062]jq
＝[0
(2m-1)q
×
(q-1-q) blkdiag{y} 0
(2m-1)q
×q] (5)
[0063]
其中，0
(2m-1)q
×
(q-1-q)
表示(2m-1)q
×
(q-1-q)维的全零矩阵，q＝0,1,2
…
,q-1，y＝[i
(q)
,b1,b2,
…
,b
2m-2
]，b1＝b2＝
…
＝b
2m-2
＝[0q×
(q-1) i
(q)
]，jq∈{0,1}
(2m-1)q
×
(2m-1)l
，blkdiag{
·
}表示生成以矩阵块为对角线的对角矩阵，q为的秩，等于分布式稀疏阵列最大可检测源信号的数量。
[0064]
进一步地，本实施方式中，步骤五中，构造toeplitz矩阵为：
[0065][0066]
其中，jq为选择矩阵，q＝0,1,2
…
,q-1。
[0067]
使用仿真实验来验证所提算法的优越性。仿真中，阵列模型选择分布式超级嵌套阵列，该分布式稀疏阵列包含2个子阵列，根据公式(5)和(6)，构造的toeplitz矩阵为：
[0068][0069]
其中，[z1]q表示z1中第q个位置的数据。分布式超级嵌套阵列的协同阵列如图3所示。
[0070]
具体仿真条件如下：源信号方位是30
°
，天线的总数等于28，快拍数为500，基本均匀间距d等于半波长，进行2000次蒙特卡罗试验。在相同的仿真条件下，分别使用空间平滑算法和本发明提出的新型算法进行doa估计，并记录两种算法的运行时长。图4为doa估计的均方根误差(rmse)随信噪比(snr)变化的仿真结果。
[0071]
图4的仿真结果表明，两种算法具有相同的估计精度，均能大大提高doa估计的精度。然而，在配备intel core i5-7500 3.40ghz处理器和8gb ram的联想计算机上，空间平滑算法的运行时间为518.24秒，本发明提出的新型算法仅花费295.77秒，节省了42.9％的运行时间。因此，当子阵列数目和天线数目特别多时，本发明提出的新型toeplitz矩阵重构算法可以大大减少运算量，提高计算速度，显著降低对计算机硬件性能的要求。仿真结果证明了所提出的新型算法的优越性。
[0072]
虽然在本文中参照了特定的实施方式来描述本发明，但是应该理解的是，这些实施例仅仅是本发明的原理和应用的示例。因此应该理解的是，可以对示例性的实施例进行许多修改，并且可以设计出其他的布置，只要不偏离所附权利要求所限定的本发明的精神和范围。应该理解的是，可以通过不同于原始权利要求所描述的方式来结合不同的从属权利要求和本文中所述的特征。还可以理解的是，结合单独实施例所描述的特征可以使用在
其他所述实施例中。