CN109521393A - 一种基于信号子空间旋转特性的波达方向估计算法 - Google Patents

Abstract

本发明涉及一种波达方向估计方法，特别是一种基于信号子空间旋转特性的波达方向估计算法，属于信号处理技术领域。本发明提出将信号子空间旋转变换思想与旋转不变子空间算法相结合，得到一种新的波达方向估计算法。本发明利用矩阵的秩亏损特性将信号子空间按行分块并进行旋转变换得到新的维度更低的信号子空间，然后结合信号子空间之间的旋转不变特性求解出波达方向角。本发明与现有技术相比，主要解决了ESPRIT算法在低信噪比和少快拍数下估计稳定性较差的问题，同时本方法利用旋转变换降低了信号子空间的维度，从而降低了算法的计算复杂度。

Description

一种基于信号子空间旋转特性的波达方向估计算法

技术领域

本发明涉及一种波达方向估计方法，特别是一种基于信号子空间旋转特性的波达方向估计算法，属于信号处理技术领域。

背景技术

波达方向估计是在空域滤波、时域谱估计的基础上发展起来的一门学科，其在雷达、被动声纳、生物医学、射电天文以及地震勘察等军事和民用领域都有着广阔的应用前景，是信号处理领域的研究热点。

信号子空间的旋转不变子空间技术(ESPRIT算法)，是利用信号子空间旋转不变特性，提取出信号的方位信息从而进行DOA估计的一种算法。该算法在高信噪比下具有良好的估计性能且避免了谱峰搜索，但在低信噪比、少快拍数时，其估计稳定性较差。

发明内容

本发明提出了一种基于信号子空间旋转特性的波达方向估计算法，主要解决了ESPRIT算法在低信噪比和少快拍数下估计稳定性较差的问题。

本发明的技术方案是：主要解决了ESPRIT算法在低信噪比和快拍数少下估计稳定性较差的问题，包括如下具体步骤：

(1)构建天线阵，包括子阵X和子阵Y

天线阵为M个阵元偶组成的均匀线阵，每个阵元偶包含两个响应特性完全相同的阵元，且阵元间距为位移矢量Δ，阵元偶之间的间距为d，将均匀线阵分为两个平移矢量为Δ的子阵X和子阵Y，子阵X由阵元偶x₁,x₂,…x_M组成，子阵Y由阵元偶y₁,y₂,…y_M组成，阵元x₁和y₁对应构成一个阵元偶，后面依次对应总共有M个阵元偶，假设P个来自空间远场窄带的独立目标信号源入射到该等距的均匀线阵上，信号的入射方向为θ₁,θ₂,…θ_P，其中，M≥P，λ为入射信号的波长。

(2)由天线阵中的天线接收机对空间信号采样分别得到子阵X和子阵Y的接收信号矢量即其输出信号矢量X(t)、Y(t)，并将其合并得到整个阵列的输出信号矢量Z(t)；

子阵X和子阵Y的输出信号矢量为：

X(t)＝[x₁(t),x₂(t),…x_M(t)]^T＝AS(t)+N_x(t)

Y(t)＝[y₁(t),y₂(t),…y_M(t)]^T＝AΦS(t)+N_y(t)

式中，表示将子阵X和子阵Y联系在一起的旋转算子，j是虚数，j²＝-1，A(t)＝[a(θ₁),a(θ₂),…,a(θ_P)]表示导向矢量，其中a(θ_i)＝exp[j(m-1)2πdsinθ_i/λ]，i＝1,2,3,…,P，a(θ_i)表示第i个信源的导向矢量，N_x和N_y分别表示子阵X和子阵Y产生的噪声矢量。

将子阵X和子阵Y的输出信号矢量合并，得到阵列总的输出信号矢量Z(t)：

(3)计算信号矢量Z(t)的自相关矩阵R_z，并对其进行特征值分解得到信号子空间；

根据阵列输出得到其自相关矩阵并进行特征值分解：

式中，R_S表示信号自相关矩阵，σ²表示噪声方差，E表示求期望，H表示矩阵的共轭转置，U_S表示特征分解后前P个较大特征值对应的特征矢量张成信号子空间，U_N表示其余的2M-P个较小特征值对应的特征矢量张成噪声子空间；特征分解后的特征值从大到小排，排在前P个的特征值等于信号和噪声方差之和，后面的2M-P个特征值等于噪声的方差，Λ_S和Λ_N分别表示P个较大特征值特征值和后面的2M-P个特征值形成的对角矩阵。

(4)对信号子空间按行进行分块处理，得到低维度的子阵信号子空间U_new；由于信号子空间U_S为2M×P维矩阵，将矩阵按行进行分块：其中，子阵U_S1为2M-P×P维矩阵，U_S2为P×P维矩阵；

由于U_S2为行满秩矩阵，则存在右逆矩阵使得I为单位阵；

将原信号子空间U_S进行变换处理：

可以看出，对原信号子空间的变换处理实际上可理解为矩阵U_S列向量的加权和，也就是将矩阵U_S的列空间span(U_S)在2M维空间进行了角度旋转。且旋转之后的矩阵满足：U_new矩阵的秩小于矩阵U_S和秩的最小值，即旋转处理后的信号子空间是原信号子空间的子集。

(5)利用子阵信号子空间之间的旋转不变特性构造旋转关系矩阵，从而求得波达方向角，由于特征矢量张成的信号子空间与信号的导向矢量张成的空间是同一空间，则经过变换处理后满足：span{U_new}＝span{a(θ₁),a(θ₂),…,a(θ_P)}，

所以，存在一个唯一的非奇异矩阵T，使得：

阵列输出由子阵X和子阵Y的输出合并得到

U_newY＝U_newXT^-1ΦT＝U_newXΨ，且当A为满秩矩阵时，Φ＝TΨT^-1；

最后，利用最小二乘法可计算出矩阵对矩阵Ψ进行特征分解，其P个特征值对应于P个信号的到达角，从而求得信号的波达方向角。

本发明的有益效果是：本发明将信号子空间旋转变换思想与旋转不变子空间算法相结合，得到一种新的波达方向算法。利用矩阵的秩亏损特性将信号子空间按行分块并进行旋转变换得到新的维度更低的信号子空间，然后利用新的信号子空间之间的旋转不变特性求解出信号的波达方向角。本发明通过旋转变换使得新的信号子空间中包含了更多的信号参数信息，从而提高了ESPRIT算法在低信噪比和快拍数少下估计稳定性，同时本方法降低了信号子空间的维度，从而降低了算法的计算复杂度。

附图说明

图1是本发明的总流程图；

图2是本发明与ESPRIT算法在不同信噪比条件下均方误差对比图；

图3是本发明与ESPRIT算法在不同采样数下均方误差对比图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明。

实施例1：如图1所示，一种基于信号子空间旋转特性的波达方向估计算法，具体步骤如下：

步骤1：利用天线接收机形成两个子阵

每间隔距离d放置一对天线接收机，每对接收机包含两个响应特性完全相同的阵元，且阵元间距为已知位移矢量Δ，一共放置M个，形成均匀线阵，将阵列分为两个平移矢量为Δ的子阵X和子阵Y，且子阵X和子阵Y分别由阵元偶的x₁,x₂,…x_M和y₁,y₂,…y_M构成。假设有P个远场窄带信号入射到均匀线阵上，且信号在传播过程中加入了零均值的高斯白噪声，其中，M≥P，λ为入射信号的波长。

步骤2：由阵列天线接收机对空间信号采样分别得到子阵X和子阵Y的接收信号矢量即输出信号矢量X(t)、Y(t)，并将输出合并得到整个阵列的输出信号矢量Z(t)。

为了描述方便，我们将阵列分为两个平移矢量为Δ的子阵X和子阵Y，且子阵X和子阵Y分别由阵元偶的x₁,x₂,…x_M和y₁,y₂,…y_M构成。

2)则两个子阵产生输出矢量为：

X(t)＝[x₁(t),x₂(t),…x_M(t)]^T＝AS(t)+N_x(t)

Y(t)＝[y₁(t),y₂(t),…y_M(t)]^T＝AΦS(t)+N_y(t)

式中，表示将子阵1和子阵2联系在一起的旋转算子，A(t)＝[a(θ₁),a(θ₂),…,a(θ_P)]表示导向矢量，其中a(θ_i)＝exp[j(m-1)2πdsinθ_i/λ]，i＝1,2,3,…,P，a(θ_i)表示第i个信源的导向矢量，N_x和N_y分别表示子阵1和子阵2产生的噪声矢量。

3)将两个子阵的输出合并，得到阵列总的输出为：

步骤3：计算信号矢量Z(t)的自相关矩阵R_z，并对其进行特征值分解得到信号子空间。

根据阵列的输出得到其自相关矩阵并进行特征值分解：

式中，R_S表示信号自相关矩阵，σ²表示噪声方差，E表示求期望，H表示矩阵的共轭转置，U_S表示特征分解后P个较大的特征值对应的特征矢量张成信号子空间，U_N表示其余的2M-P个较小特征值对应的特征矢量张成噪声子空间；Λ_S和Λ_N分别表示特征值形成的对角矩阵。

步骤4：对信号子空间按行进行分块处理，得到低维度的子阵信号子空间U_new。

1)由于信号子空间U_S为2M×P维矩阵，可将矩阵按行进行分块：其中，子阵U_S1为(2M-P)×P维矩阵，U_S2为P×P维矩阵。

2)由于U_S2为行满秩矩阵，即存在右逆矩阵满足公式

3)将原信号子空间进行变换处理：

4)可以看出，对原信号子空间的变换处理实际上可理解为矩阵U_S列向量的加权和，也就是将矩阵U_S的列空间span(U_S)在2M维空间进行了角度旋转。且旋转之后的矩阵满足：

∴即旋转处理后的信号子空间是原信号子空间的子集。

5)由于特征矢量张成的信号子空间与信号的导向矢量张成的空间是同一空间，则经过变换处理后满足：span{U_new}＝span{a(θ₁),a(θ₂),…,a(θ_P)}

所以，存在一个唯一的、非奇异矩阵T，使得：

6)∵阵列输出有两个子阵的输出合并得到

∴

∴U_newY＝U_newXT^-1ΦT＝U_newXΨ，且当A为满秩矩阵时，Φ＝TΨT^-1。

7)最后，利用最小二乘法可计算出矩阵对矩

阵Ψ进行特征分解，其P个特征值对应于P个信号的到达角，从而求得信号的波达方向角。

本发明利用矩阵的秩亏损特性将信号子空间按行分块并进行旋转变换得到新的信号子空间，然后利用新的信号子空间之间的旋转不变特性求解出信号的波达方向角。该算法通过旋转变换，使得新的信号子空间中包含了更多的信号参数信息，提高了ESPRIT算法在低信噪比和快拍数少下估计稳定性，同时本方法利用旋转变换降低了信号子空间的维度，从而降低了算法的计算复杂度。

实施例2：按照实施例1中的方法来进行计算，其中考虑8个全向阵元组成的均匀线阵，子阵阵元数为7个，将前7个阵元输出表示为子阵1，后7个阵元输出表示为子阵2，阵元间距为0.5，空域角度搜索范围为[-90° 90°]。

均方误差计算式为：

式中，I表示蒙特卡洛实验次数，表示第i次试验的波达方向角，θ_p表示信号真实的波达方向角。

仿真1：假设三个独立信源分别以角度[-60° 0° 30°]入射到该等距均匀线阵上，采样点数为100，使信噪比由-10dB增加到10dB，每个信噪比下都进行50次独立的波达方向估计试验，利用本发明算法与ESPRIT算法进行对比，分别计算不同信噪比下两种算法的均方误差，试验结果如图2所示，其中：图2中横坐标表示信噪比，纵坐标表示均方误差。

仿真2：在仿真1的条件下，设置信噪比为-10dB，使采样数由100增加到200，每个采样点数下进行50次独立的波达方向估计试验，将本发明算法与ESPRIT算法进行对比，分别计算不同采样数下两种算法的均方误差，试验结果如图3所示，其中：图3中横坐标表示采样数，纵坐标表示均方误差。

从图2、图3中可以看出，随着信噪比和采样数的增加两种算法波达方向角的估计均方误差出现波动，但是在低信噪比、少快拍数情况下，本发明的均方误差波动明显小于ESPRIT算法，而且本发明算法的均方误差值也小于ESPRIT算法，说明了本方法的稳定性和有效性优势。

综上所述，本发明算法在信噪比低和少快拍数下估计性能优良、具有很好的稳定性和有效性。

上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

Claims (3)

1.一种基于信号子空间旋转特性的波达方向估计算法，其特征在于：具体步骤如下：

(1)构建天线阵，包括子阵X和子阵Y

天线阵为M个阵元偶组成的均匀线阵，每个阵元偶包含两个响应特性完全相同的阵元，且阵元间距为位移矢量Δ，阵元偶之间的间距为d，将均匀线阵分为两个平移矢量为Δ的子阵X和子阵Y，子阵X由阵元偶x₁,x₂,…x_M组成，子阵Y由阵元偶y₁,y₂,…y_M组成，假设P个来自空间远场窄带的独立目标信号源入射到该等距的均匀线阵上，信号的入射方向为θ₁,θ₂,…θ_P，其中，M≥P，λ为入射信号的波长；

(2)由天线阵中的天线接收机对空间信号采样分别得到子阵X和子阵Y的接收信号矢量即其输出信号矢量X(t)、Y(t)，并将其合并得到整个阵列的输出信号矢量Z(t)；

(3)计算信号矢量Z(t)的自相关矩阵R_z，并对其进行特征值分解得到信号子空间；

(4)对信号子空间按行进行分块处理，得到低维度的子阵信号子空间U_new；

(5)利用子阵信号子空间之间的旋转不变特性构造旋转关系矩阵，从而求得波达方向角。

2.根据权利要求1所述的基于信号子空间旋转特性的波达方向估计算法，其特征在于：所述步骤(2)和步骤(3)的具体过程如下：

子阵X和子阵Y的输出信号矢量为：

X(t)＝[x₁(t),x₂(t),…x_M(t)]^T＝AS(t)+N_x(t)

Y(t)＝[y₁(t),y₂(t),…y_M(t)]^T＝AΦS(t)+N_y(t)

式中，表示将子阵X和子阵Y联系在一起的旋转算子，j是虚数，A(t)＝[a(θ₁),a(θ₂),…,a(θ_P)]表示导向矢量，其中a(θ_i)＝exp[j(m-1)2πd sinθ_i/λ]，i＝1,2,3,…,P，a(θ_i)表示第i个信源的导向矢量，N_x和N_y分别表示子阵X和子阵Y产生的噪声矢量；

将子阵X和子阵Y的输出信号矢量合并，得到阵列总的输出信号矢量Z(t)：

根据阵列输出得到其自相关矩阵并进行特征值分解：

式中，R_S表示信号自相关矩阵，σ²表示噪声方差，E表示求期望，H表示矩阵的共轭转置，U_S表示特征分解后前P个较大特征值对应的特征矢量张成信号子空间，U_N表示其余的2M-P个较小特征值对应的特征矢量张成噪声子空间；Λ_S和Λ_N分别表示特征值形成的对角矩阵。

3.根据权利要求1所述的基于信号子空间旋转特性的波达方向估计算法，其中步骤(4)和步骤(5)的具体过程如下：

a)信号子空间的分块处理

由于信号子空间U_S为2M×P维矩阵，将矩阵按行进行分块：其中，子阵U_S1为(2M-P)×P维矩阵，U_S2为P×P维矩阵；

由于U_S2为行满秩矩阵，则存在右逆矩阵使得I为单位阵；

将原信号子空间U_S进行变换处理：

处理后的矩阵满足：U_new矩阵的秩小于矩阵U_S和秩的最小值，即旋转处理后的信号子空间是原信号子空间的子集；

b)利用子阵的旋转不变关系求解信号的波达方向角

由于特征矢量张成的信号子空间与信号的导向矢量张成的空间是同一空间，则经过变换处理后满足：span{U_new}＝span{a(θ₁),a(θ₂),…,a(θ_P)}，

所以，存在一个唯一的非奇异矩阵T，使得：

阵列输出由子阵X和子阵Y的输出合并得到

U_newY＝U_newXT^-1ΦT＝U_newXΨ，且当A为满秩矩阵时，Φ＝TΨT^-1；

最后，利用最小二乘法可计算出矩阵对矩阵Ψ进行特征分解，其P个特征值对应于P个信号的到达角，从而求得信号的波达方向角。