CN105388462B - 一种基于互质阵列的虚拟波束形成方法 - Google Patents

Abstract

本发明公开了一种基于互质阵列的虚拟波束形成方法，先布设多个阵元形成互质阵列，对不同空间方向上多个非相关信号源的入射信号进行多快拍采样；再对不同相关间隔的相关函数值进行时间平均估计；再基于相关间隔与互质阵列差协同阵虚拟阵元位置的对应关系，提取具有连续相关间隔的相关函数值，构造差协同阵中均匀密布虚拟阵元的等效单快拍信号；最后根据均匀密布虚拟阵元位置分布设计权向量，对等效的单快拍信号进行加权叠加得到输出波束。本发明利用了互质阵列的差协同阵中大部分虚拟阵元均匀密布的特点，解决了现有稀疏子阵列波束联乘方法存在的旁瓣干扰、不适用于多信号源的问题，并且有效地提高了波束信噪比，增强了抗噪声能力。

Description

一种基于互质阵列的虚拟波束形成方法

技术领域

本发明属于阵列波束形成技术领域，涉及一种基于互质阵列的虚拟波束形成方法。

背景技术

互质阵列是在2010年提出的一种新型的非均匀稀疏阵列，由两个阵元间距呈现互质关系的均匀稀疏子阵列组成，相比于最小冗余矩阵、嵌套阵列等传统非均匀稀疏阵列，互质阵列具有阵元布局设计简单、阵元互耦影响微小等突出优点。

阵列信号处理在雷达、声呐、通信、地震勘探等领域获得了广泛应用，它的两个主要研究方向为波达方向估计和波束形成。互质阵列最早应用于空间信号源的波达方向估计，利用差协同阵等效原理扩张相关矩阵的维数，提高可分辨信号源的数目，实现最大可分辨信号源数目突破阵元数目的限制。

互质阵列也逐渐被应用于波束形成方面，现有的方法概括为两步：第一步，两个均匀稀疏子阵列的接收信号分别进行加权叠加，通过权系数沿空间方向的遍历，形成两个独立的带有栅瓣的子波束；第二步，两个子波束进行联乘融合，抑制栅瓣同时保留在信号源方向的主瓣，达到提取期望信号的目的。然而，该方法在子波束联乘融合时部分栅瓣被保留退化为旁瓣，造成明显的旁瓣干扰，此外，空间中存在多个信号源时，子波束中信号源的栅瓣会混叠入主瓣，导致子波束主瓣扭曲或衰退，联乘融合造成严重的栅瓣残余和主瓣退化，输出波束杂乱无章无法分辨多个信号源方向。因此，现有基于互质阵列的波束形成方法存在旁瓣干扰、不适用于多信号源的明显缺陷。

发明内容

本发明的目的就在于提供一种解决上述问题，可以明显克服旁瓣干扰、能够轻易分辨出多信号源的基于互质阵列的虚拟波束形成方法。

为了实现上述目的，本发明采用的技术方案是这样的：

一种基于互质阵列的虚拟波束形成方法，包括以下步骤：

1)布设多个阵元形成互质阵列，对不同空间方向上多个非相关信号源的入射信号进行多快拍采样；

2)利用互质阵列采集的多快拍信号，对互质阵列的相关矩阵中，不同相关间隔的相关函数值进行时间平均估计；

3)基于相关间隔与互质阵列差协同阵虚拟阵元位置的对应关系，提取具有连续相关间隔的相关函数值，构造差协同阵中均匀密布虚拟阵元的等效单快拍信号；

4)根据均匀密布虚拟阵元位置分布设计权向量，对等效的单快拍信号进行加权叠加得到输出波束。

作为优选：步骤1)具体方式为：将N+2M－1个阵元摆放为包含两个均匀稀疏子阵列的互质阵列，阵元位置为：

x＝{Mnd₀,0≤n≤N-1}U{Nmd₀,1≤m≤2M-1} (1)

式(1)中，d₀为半波长的单位间隔，当工作频率为f₀时d₀＝λ/2＝c/2f₀；在互质阵列的差协同阵中，虚拟阵元的位置由式(1)中阵元位置的差值决定，中间2MN+2M-1个虚拟阵元以单位间隔d₀为阵元间距均匀分布于-(MN+M-1)d₀到(MN+M-1)d₀之间，即中间2MN+2M-1个虚拟阵元虚拟阵元呈现均匀密布特性，而边缘两侧小部分的虚拟阵元呈现非均匀稀疏分布特性，如图2所示。

当功率分别为L,的L个互不相关的信号源沿各自对应方向θ₁,θ₂,L,θ_L入射至互质阵列，则在第k个快拍采样时刻的接收信号矢量表示为：

y(k)＝[y₁(k),y₂(k),L,y_N+2M-1(k)]

(2)

＝A×s(k)+ε(k)

式(2)中，s(k)＝[s₁(k),s₂(k),L,s_L(k)]^T为源信号矢量，上标T表示矩阵转置；ε(k)为空域和时域的白噪声矢量；A是互质阵列操纵矩阵，A中第(i,l)个元素表示为：

式(3)中，x_i为互质阵列中第i个阵元的位置，θ_j是第j个信号源的入射方向，k₀＝2πf₀/c为波数，c为自由空间光速；

若互质阵列共进行了K次快拍采样，形成的信号矩阵表示为：

Y＝[y(1),y(2),L,y(K)] (4)。

作为优选：步骤2)具体为：

假设互质阵列中第q个阵元与第i个阵元的间距为半波长单位间隔的p倍，即x_q-x_i＝pd₀，同时也是差协同阵中某个虚拟阵元的位置，则相关矩阵R的时间平均估计中，与虚拟阵元位置pd₀对应的相关间隔为p的相关函数值计算为：

式(5)中，*表示共轭。

理论上，互质阵列的相关矩阵计算为：

式(9)中，为源信号相关矩阵，为噪声功率，I为单位矩阵。上标H表示转置共轭操作，E[·]为统计平均操作符，相关矩阵R中相关间隔为p的相关函数值表示为：

式(10)中，pd₀为互质阵列两个阵元的间隔，即差协同阵中虚拟阵元位置，与相关间隔p一一对应，且相关间隔等于虚拟阵元位置除以半波长单位间隔d₀。

式(9)中相关矩阵的统计平均计算需要互质阵列信号的所有样本，实际情况下无法实现，因此，本发明采用式(4)中K次快拍的有限样本计算相关矩阵R的时间平均估计假设互质阵列中第q个阵元与第i个阵元的间距为半波长单位间隔的p倍，即x_q-x_i＝pd₀，则对应中相关间隔为p的相关函数值(即式(10)中r(p)的时间平均估计)计算为：

作为优选：步骤3)具体为，基于相关间隔与虚拟阵元位置的一一对应关系，从中提取相关间隔从-(MN-M+1)至MN-M+1连续变化的相关函数值，将它们顺序排列则构造出差协同阵中分布于-(MN+M-1)d₀到(MN+M-1)d₀之间的2MN+2M-1个均匀密布虚拟阵元的等效单快拍信号矢量，即：

根据式(10)中相关函数值的表达式，式(6)中的信号矢量等效于互质阵列差协同阵中2MN+2M-1个均匀密布虚拟阵元，对空间方向θ₁,θ₂,L,θ_L上L个相关信号源入射信号的单快拍采样，且L个信号源的入射信号幅度分别为L,

作为优选：步骤4)具体为，

根据-(MN+M-1)d₀到(MN+M-1)d₀的均匀密布虚拟阵元的位置，虚拟波束形成的权向量设计为：

W＝[e^{j(MN+M-1)πsinθ},L,e^j(1)πsinθ,1,e^-j(1)πsinθ,L e^{-j(MN+M-1)πsinθ}] (7)

则针对式(6)中虚拟阵元等效单快拍信号矢量，虚拟波束形成表示为：

其中θ的取值需以间隔Δθ渐进递增，遍历-90°至90°的空间方向区间后得到最终的输出波束。

与现有技术相比，本发明的优点在于：本发明提出的基于互质阵列的差协同阵虚拟波束形成方法，充分利用了互质阵列的差协同阵中大部分虚拟阵元均匀密布的特点，将对应的连续相关间隔的相关函数值构造为等效的单快拍信号进行虚拟波束形成。

相比于现有的稀疏子阵列波束联乘方法，该方法解决了旁瓣干扰、不适用于多信号源的问题，由波束主瓣指向最多可分辨MN+M－1个信号源方向，此外，本发明在估计相关函数值时使用了互质阵列的多快拍信号，相比于现有利用单快拍信号的方法，有效地提高了波束信噪比，增强了抗噪声能力。

附图说明

图1为本发明虚拟波束形成方法的处理流程图；

图2为本发明中互质阵列及其差协同阵的布局示意图；

图3为本发明实施例1中互质阵列及其差协同阵的布局示意图；

图4为空间有单个信号源时，本发明虚拟波束形成示意图；

图5为空间有单个信号源时，互质阵列中子阵列波束波形图；

图6为图5中两个子阵列波束联乘后波形图；

图7为空间有三个信号源时，本发明虚拟波束形成示意图；

图8为空间有三个信号源时，互质阵列中子阵列波束波形图；

图9为图8中两个子阵列波束联乘后波形图；

图10为空间有七个信号源时，本发明虚拟波束形成示意图；

图11为空间有七个信号源时，互质阵列中子阵列波束波形图；

图12为图11中两个子阵列波束联乘后波形图。

具体实施方式

下面将结合附图对本发明作进一步说明。

实施例1：参见图1到图12，一种基于互质阵列的虚拟波束形成方法，包括以下步骤：

1)布设多个阵元形成互质阵列，对不同空间方向上多个非相关信号源的入射信号进行多快拍采样；具体方式为：将N+2M－1个阵元摆放为包含两个均匀稀疏子阵列的互质阵列，阵元位置为：

x＝{Mnd₀,0≤n≤N-1}U{Nmd₀,1≤m≤2M-1} (1)

式(1)中，d₀为半波长的单位间隔，当工作频率为f₀时d₀＝λ/2＝c/2f₀；本实施例中，我们设置M＝2、N＝3，半波长单位间隔d₀＝1。

当功率分别为L,的L个互不相关的信号源沿各自对应方向θ₁,θ₂,L,θ_L入射至互质阵列，则在第k个快拍采样时刻的接收信号矢量表示为：

y(k)＝[y₁(k),y₂(k),L,y_N+2M-1(k)]

(2)

＝A×s(k)+ε(k)

式(2)中，s(k)＝[s₁(k),s₂(k),L,s_L(k)]^T为源信号矢量，上标T表示矩阵转置；ε(k)为空域和时域的白噪声矢量；A是互质阵列操纵矩阵，A中第(i,l)个元素表示为：

式(3)中，x_i为互质阵列中第i个阵元的位置，θ_j是第j个信号源的波达方向，k₀＝2πf₀/c为波数，c为自由空间光速；

若互质阵列共进行了K次快拍采样，形成的信号矩阵表示为：

Y＝[y(1),y(2),L,y(K)] (4)；

本实施例中，我们预设互质阵列的每个阵元均进行了1000次快拍采样；

2)利用互质阵列采集的多快拍信号，对互质阵列的相关矩阵中，不同相关间隔的相关函数值进行时间平均估计，具体为：

假设互质阵列中第q个阵元与第i个阵元的间距为半波长单位间隔的p倍，即x_q-x_i＝pd₀，其反映在差协同阵中某个虚拟阵元的位置，则相关矩阵R的时间平均估计中，与虚拟阵元位置pd₀对应的相关间隔为p的相关函数值计算为：

理论上，互质阵列的相关矩阵计算为：

其中，为源信号相关矩阵，为噪声功率，I为单位矩阵。上标H表示转置共轭操作，E[·]为统计平均操作符。相关矩阵R中相关间隔为p的相关函数值表示为：

其中，pd₀为互质阵列两个阵元的间隔，即差协同阵中虚拟阵元位置，与相关间隔p一一对应，且相关间隔等于虚拟阵元位置除以半波长单位间隔d₀。

式(9)中相关矩阵的统计平均计算需要互质阵列信号的所有样本，实际情况下无法实现，因此，本发明采用式(4)中K次快拍的有限样本计算相关矩阵R的时间平均估计假设互质阵列中第q个阵元与第i个阵元的间距为半波长单位间隔的p倍，即x_q-x_i＝pd₀，则中相关间隔为p的相关函数值(即式(10)中r(p)的时间平均估计)计算为：

3)基于相关间隔与互质阵列差协同阵虚拟阵元位置的对应关系，提取具有连续相关间隔的相关函数值，构造差协同阵中均匀密布虚拟阵元的等效单快拍信号；

具体为，从中提取相关间隔从-(MN-M+1)至MN-M+1连续变化的相关函数值，将它们顺序排列则构造出差协同阵中分布于-(MN+M-1)d₀到(MN+M-1)d₀之间的2MN+2M-1个均匀密布虚拟阵元的等效单快拍信号矢量，即：

4)根据均匀密布虚拟阵元位置分布设计权向量，对等效的单快拍信号进行加权叠加得到输出波束；

具体为，根据均匀密布虚拟阵元的位置分布，虚拟波束形成的权向量设计为：

W＝[e^{j(MN+M-1)πsinθ},L,e^j(1)πsinθ,1,e^-j(1)πsinθ,L e^{-j(MN+M-1)πsinθ}] (7)

则虚拟波束形成表示为：

其中θ的取值需以间隔Δθ渐进递增，遍历-90°至90°的空间方向区间后得到最终的输出波束。

为了更清晰的说明上述流程，我们采用MATLAB进行仿真。

步骤1)中，仿真采用M＝2、N＝3的互质阵列，设置半波长单位间隔d₀＝1，则互质阵列的布局以及对应差协同阵的布局如图3所示，差协同阵中除去边缘处于±9位置上的两个虚拟阵元，中间的2MN+2M-1＝15个虚拟阵元以单位间隔1均匀密布于-(MN+M-1)＝-7至MN+M-1＝7的位置区间。互质阵列的每个阵元均进行了1000次快拍采样，分别设置单个信号源和多个不同方向的非相关信号源进行仿真，并加入加性高斯白噪声，信噪比设置为10dB。具体可参见图2、图3。

步骤2)中，基于互质阵列1000次快拍采样的数据，相关矩阵的时间平均估计中相关间隔为p的相关函数值计算为：

其中，互质阵列中第q个阵元与第i个阵元的位置差为x_q-x_i＝p。由于相关间隔与差协同阵虚拟阵元的位置(互质阵列两个阵元位置差)一一对应，因此，对于仿真采用的M＝2、N＝3的互质阵列，相关矩阵中相关间隔p的取值包括差协同阵中两个边缘虚拟阵元位置(±9)以及中间15个均匀密布虚拟阵元位置(﹣7至7，步进为1)。

步骤3)中，从步骤2)获得的相关矩阵估计中提取相关间隔从﹣7到7的相关函数值，依次排列构造出差协同阵中中间15个均匀密布虚拟阵元的等效单快拍信号，即：

步骤4)中，虚拟波束形成的权向量设计为：

W＝[e^j7πsinθ,L,e^jπsinθ,1,e^-jπsinθ,L,e^-j7πsinθ] (13)

则针对式(12)中均匀密布虚拟阵元的等效快拍信号矢量，虚拟波束形成表示为：

其中θ的取值从﹣90°开始以间隔Δθ＝1°渐进递增直至90°，遍历﹣90°至90°的空间方向区间后得到最终的输出波束。

图4给出了空间中存在单个信号源时利用本发明形成的虚拟波束，信号源的空间方向为0°，即以0°入射到互质阵列，经本发明方法仿真出来的虚拟波束，其主瓣指向为0°与信号源的空间方向一致，说明了本发明对单个信号源波束形成的有效性。

作为对比，现有的互质阵列子阵列波束联乘方法的仿真结果也一并给出，图5显示了两个稀疏子阵列的波束，明显观察到栅瓣的存在，图6的子阵列波束联乘结果中存在子波束栅瓣联乘融合造成的旁瓣干扰。比较而言，图4中本发明形成的虚拟波束未见明显的旁瓣干扰，证明了本发明在旁瓣干扰方面的优势。

当空间中存在多个非相干的信号源时，首先仿真设置三个信号源，信号源的空间方向分别为-30°、0°和45°，图7给出了本发明的虚拟波束形成结果，图8和图9是现有方法的子阵列波束及其联乘融合结果。对比图7和图9可知，本发明能够形成主瓣指向信号源方向的清晰波束，证明了对多个信号源波束形成的适用性，而现有子阵列波束联乘方法不再适用于多个信号源的情况，这是由于图8子波束中存在由栅瓣与主瓣混叠引起的复杂栅瓣干扰和主瓣扭曲退化而造成的。

图10至图12给出了空间中存在七个非相干信号源时的波束形成结果，首先仿真设置七个信号源，信号源的空间方向分别为﹣60°、﹣36°、﹣18°、0°、18°、36°、60°，其中图10为本发明形成的虚拟波束，可以清晰观察到七个主瓣，且主瓣指向与信号源方向一致，验证了本发明利用包含N+2M-1＝6个阵元的互质阵列最多可形成MN+M-1＝7个主瓣，即最多可分辨MN+M-1＝7个信号源的能力。而图11和图12中无法分辨指向为七个信号源方向的七个主瓣，再次证明了现有子阵列波束联乘方法对多个信号源波束形成的不适用性。

Claims (4)

1.一种基于互质阵列的虚拟波束形成方法，其特征在于：包括以下步骤：

1)布设多个阵元形成互质阵列，对不同空间方向上多个非相关信号源的入射信号进行多快拍采样；

2)利用互质阵列采集的多快拍信号，对互质阵列的相关矩阵中，不同相关间隔的相关函数值进行时间平均估计；

3)基于相关间隔与互质阵列差协同阵虚拟阵元位置的对应关系，提取具有连续相关间隔的相关函数值，构造差协同阵中均匀密布虚拟阵元的等效单快拍信号；

4)根据均匀密布虚拟阵元位置分布设计权向量，对等效的单快拍信号进行加权叠加得到输出波束；

其中，步骤1)具体方式为：将N+2M－1个阵元摆放为包含两个均匀稀疏子阵列的互质阵列，阵元位置为：

x＝{Mnd₀,0≤n≤N-1}∪{Nmd₀,1≤m≤2M-1} (1)

式(1)中，d₀为半波长的单位间隔，当工作频率为f₀时d₀＝λ/2＝c/2f₀；

当功率分别为的L个互不相关的信号源沿各自对应方向θ₁,θ₂,…,θ_L入射至互质阵列，则在第k个快拍采样时刻的接收信号矢量表示为：

$\begin{array}{c}y\left(k\right)=\[y_1\left(k\right),y_2\left(k\right),...,y_{N+2M-1}\left(k\right)\]\\ =A\×s\left(k\right)+\mathrm{\ϵ}\left(k\right)\end{array}---\left(2\right)$

式(2)中，s(k)＝[s₁(k),s₂(k),…,s_L(k)]^T为源信号矢量，上标T表示矩阵转置；ε(k)为空域和时域的白噪声矢量；A是互质阵列操纵矩阵，A中第(i,l)个元素表示为：

${\[A\]}_{i,l}=e^{{\mathrm{jk}}_0{x}_{i}sin\left({\mathrm{\θ}}_l\right)},i=1,2,...,N+2M-1;l=1,2,...,L---\left(3\right)$

式(3)中，x_i为互质阵列中第i个阵元的位置，θ_j是第j个信号源的入射方向，k₀＝2πf₀/c为波数，c为自由空间光速；

若互质阵列共进行了K次快拍采样，形成的信号矩阵表示为：

Y＝[y(1),y(2),…,y(K)] (4)。

2.根据权利要求1所述的一种基于互质阵列的虚拟波束形成方法，其特征在于：步骤2)具体为，

假设互质阵列中第q个阵元与第i个阵元的间距为半波长单位间隔的p倍，即x_q-x_i＝pd₀，则相关矩阵R的时间平均估计中，相关间隔为p的相关函数值计算为：

$\hat{r}\left(p\right)=\frac{1}{K}\underset{k=1}{\overset{K}{\Σ}}{y}_i\left(k\right)y_q^{*}\left(k\right)---\left(5\right).$

3.根据权利要求2所述的一种基于互质阵列的虚拟波束形成方法，其特征在于：步骤3)具体为，从中提取相关间隔从-(MN-M+1)至MN-M+1连续变化的相关函数值，将它们顺序排列则构造出差协同阵中2MN+2M-1个均匀密布虚拟阵元的等效单快拍信号矢量，即：

$z=\[\hat{r}(-MN-M+1),...,\hat{r}(-1),\hat{r}\left(0\right),\hat{r}\left(1\right),...,\hat{r}(MN+M-1)\]---\left(6\right).$

4.根据权利要求3所述的一种基于互质阵列的虚拟波束形成方法，其特征在于：步骤4)具体为，

根据均匀密布虚拟阵元的位置分布，虚拟波束形成的权向量设计为：

W＝[e^{j(MN+M-1)πsinθ},…,e^j(1)πsinθ,1,e^-j(1)πsinθ,…e^{-j(MN+M-1)πsinθ}] (7)

则虚拟波束形成表示为：

$H\left(\mathrm{\θ}\right)=W\×z^T=\underset{p=-MN-M+1}{\overset{MN+M-1}{\Σ}}\hat{r}\left(p\right)e^{-jp\mathrm{\π}sin\mathrm{\θ}}---\left(8\right);$

其中θ的取值需以间隔Δθ渐进递增，遍历-90°至90°的空间方向区间后得到最终的输出波束。