CN107894581A

CN107894581A - 一种宽带阵列波达方向估计方法

Info

Publication number: CN107894581A
Application number: CN201711135852.1A
Authority: CN
Inventors: 沈明威; 张琪; 王冠; 吉雨; 王晓冬; 姚旭
Original assignee: Hohai University HHU
Current assignee: Hohai University HHU
Priority date: 2017-11-16
Filing date: 2017-11-16
Publication date: 2018-04-10

Abstract

本发明公开了一种宽带阵列波达方向（Direction of Arrival，DOA）估计方法，在原有投影子空间正交性测试算法（Test of Orthogonality of Projected Subspaces,TOPS）的基础上进行改进，将信号子空间投影到信号增强区，并利用点积向量的长度代替原来使用奇异值分解的方式衡量信号子空间与噪声子空间的正交性，从而计算出信号角度的估计值。仿真结果表明：与传统TOPS算法相比，本发明解决了原算法中的伪峰问题，提高了角度估计精度，同时有效降低了算法的计算量，提高了系统的运算效率。

Description

一种宽带阵列波达方向估计方法

技术领域

本发明涉及阵列信号角度估计领域，尤其涉及一种宽带阵列波达方向估计方法，具体为一种基于改进投影子空间正交性测试算法(Test of Orthogonality of ProjectedSubspaces,TOPS)的宽带阵列波达方向估计(Direction of Arrival,DOA)方法。

背景技术

DOA估计技术是现代阵列信号处理中的研究热点之一，主要用来估计空间信号的参数，如信号入射方向，频率等。宽带阵列信号的DOA估计算法主要有最大似然估计类方法和子空间类方法。最大似然类估计方法虽然能得到最优解，但是运算复杂，难以实现工程运用。子空间类方法主要有非相干信号子空间ISM算法和相干信号子空间CSM算法。ISM算法将宽带信号分解为若干窄带信号，分别计算各个频率点上的信号协方差矩阵，然后综合每个频点的处理结果得到整个频带范围的DOA估计。CSM算法通过选取聚焦频率，构造聚焦矩阵，将各个不同频点上的数据聚焦到一个已知频率上。CSM算法比较灵活，但是性能取决于聚焦矩阵，因此不能保证算法的稳健性。

TOPS算法是最近新提出来的一种宽带DOA方法。该算法首先将所有频点的协方差矩阵进行特征值分解，然后通过对角变换矩阵的作用，在不改变方位信息的前提下，将各频点的阵列导引矢量均用参考频率点的阵列导引矢量表示出来，最后利用各频点信号子空间和噪声子空间之间的正交性进行宽带信号DOA估计。TOPS算法不需要初始的聚焦角度，因而避免了相干信号子空间算法中由角度预估计导致的误差。但是，算法性能受参考频率点的信号子空间的估计误差影响，利用了投影矩阵修正也无可避免出现伪峰现象。

发明内容

本发明所要解决的技术问题是针对背景技术中所涉及到的缺陷，提供一种宽带阵列波达方向估计方法。

本发明为解决上述技术问题采用以下技术方案：

一种宽带阵列波达方向估计方法，包括如下步骤：

步骤1)，令接收机天线阵列为均匀直线阵列，阵元数为M，第m个阵元第n个时刻的接收信号其中，1≤n≤N，1≤m≤M，P为宽带信源个数，s_p(n)为第p个信源在第一个阵元处第n个时刻的接收信号，为第m个阵元相对于第一个阵元接收信号的时间延迟，d为阵元间距，c为光速，θ_p为第p个信源的DOA角度，n_m(n)为第m个阵元第n个时刻接收的噪声，N为阵元接收信号时域采样点数；

步骤2)，将阵列接收信号分Q段分别进行处理，Q为预先设定的小于等于N的自然数，且为整数，第m个阵元在子段q内第k个时刻的接收信号X_{m_q}(k)为:

对子段q内第m个阵元的接收信号进行FFT变换，得到其在第j个频点的频域X_{m_q}(ω_j)：

其中，ω_j为阵元接收信号FFT变换后第j个频点处的信号频率，ω_j∈[ω_L,ω_H]，j＝1，2...，J为子段q内FFT变换的频点数，ω_L、ω_H分别为子段q内第m个阵元接收信号的最低频率与最高频率，且ω_L＝ω₁，ω_H＝ω_J；S_{p_q}(ω_j)为子段q内第p个信源在第一个阵元处的接收信号FFT变换后第j个频点处的频域表达式，N_{m_p}(ω_j)为子段q内第m个阵元接收的噪声FFT变换后第j个频点处的频域表达式；

子段q内阵列接收信号的频域表达式为x_q(ω_j)＝A_q(ω_j,θ)s_q(ω_j)+n_q(ω_j)；

其中，x_q(ω_j)＝[X_{1_q}(ω_j),X_{2_q}(ω_j)，...，X_{M_q}(ω_j)]^T，s_q(ω_j)＝[S_{1_q}(ω_j)，S_{2_q}(ω_j)，...，S_{P_q}(ω_j)]^T，n_q(ω_j)＝[N_{1_q}(ω_j)，N_{2_q}(ω_j)，...，N_{M_q}(ω_j)]^T，θ＝[θ₁，θ₂，…，θ_P]^T，A_q(ω_j，θ)是M×P维关于波达角度和频率的矩阵，在第j个频点的表达式为

步骤3)，计算信号子空间和噪声子空间：

第j个频点处阵列接收信号的协方差矩阵为j＝1，2，...，J；

对R(ω_j)进行特征值分解，得到特征向量其中按其对应的特征值从大到小排列，对应的信号子空间和噪声子空间分别为

步骤4)，定义对角变换矩阵，计算任意频点对应的信号子空间：

步骤4.1)，在J个频点中等间距选取Z＝9个频点进行计算，并选取其中一个频点为参考频点，令参考频点频率为ω₀，定义对角变换矩阵z＝1，2，...，Z-1，ω_z为选取的第z个频点频率，频率差Δω_z＝ω_z-ω₀，φ为待估计的波达角度；

步骤4.2)，利用对角矩阵Φ(Δω_z，φ)将参考频点ω₀对应的信号子空间转化为任意频点对应的信号子空间，则第z个频点对应的信号子空间U_z(θ)为：

U_z(θ)＝Φ(Δω_z，φ)F₀

式中，z＝1，2，...，Z-1，F₀为参考频点ω₀对应的信号子空间；

步骤5)，计算信号子空间与噪声子空间正交性测试矩阵：

步骤5.1)，定义第z个频点处的阵列导引矢量z＝1，2，...，Z-1，其正交投影矩阵I为单位矩阵；

步骤5.2)，将信号子空间投影到信号增强区得到为第z个频点对应的投影后的信号子空间；

步骤5.3)，定义矩阵D^Im(φ)为信号子空间与噪声子空间正交性测试矩阵，W_z为第z个频点对应的噪声子空间；

步骤6)，在阵列视角范围内进行空间角度搜索，波达角度的估计值位于下式谱峰位置：

作为本发明一种宽带阵列波达方向估计方法进一步的优化方案，所述步骤2)中，

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

1.通过将信号子空间投影到信号增强区，并利用点积代替原来特征值分解的方式衡量信号子空间与噪声子空间的正交性，解决了原来TOPS算法虚假谱峰的问题。

2.提高了角度估计的精度，能比较精确的估计两个信源相互靠近的宽带信号DOA信息。

3.有效降低了算法的计算量，提高系统运算效率，利于工程实施。

附图说明

图1为本发明信号处理流程图；

图2为阵列信号的几何构型；

图3(a)为信噪比为15dB时，TOPS算法与改进TOPS算法的性能对比图；

图3(b)为信噪比为25dB时,TOPS算法与改进TOPS算法的性能对比图；

图4为信源靠近时，TOPS算法与改进TOPS算法的性能对比图；

图5为DOA估计均方根误差与信噪比关系；

图6(a)为频点数设为7时，TOPS算法与改进TOPS算法计算复杂度对比图；

图6(b)为阵元数设为8时，TOPS算法与改进TOPS算法计算复杂度对比图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

下图1是本发明的信号处理流程图。

1、接收阵列的几何构型如图2所示。假设接收机天线阵列为均匀直线阵列，阵元数为M，相邻阵元间的距离为d，天线阵基线在x轴，法线方向与阵列的方向的角度为90°，信号源为远场宽带信号，信源数为P，第p个信源的波达角度为θ_p。以第一个阵元为参考阵元，信号到达第m个阵元的时延为：τ_m＝(m-1)dsinθ/c，m＝1，...，M，c表示光速。

第m个阵元第n个时刻的接收信号为：

s_p(n)为第p个信源在第一个阵元处第n个时刻的接收信号，n_m(n)为第m个阵元第n个时刻接收的噪声，N为阵元接收信号时域采样点数。

2、阵列接收信号分Q段分别进行处理，可以取则第m个阵元在子段q内第k个时刻的接收信号X_{m_q}(k)为：

对子段q内第m个阵元的接收信号进行FFT变换，可以得到其在第j个频点的频域表达式：

其中ω_j为阵元接收信号FFT变换后第j个频点处的信号频率，ω_j∈[ω_L,ω_H]，j＝1,2...,J，为子段q内FFT变换的频点数，ω_L,ω_H分别为子段q内第m个阵元接收信号的最低频率与最高频率，且ω_L＝ω₁，ω_H＝ω_J。S_{p_q}(ω_j)为子段q内第p个信源在第一个阵元处的接收信号FFT变换后第j个频点处的频域表达式，N_{m_p}(ω_j)为子段q内第m个阵元接收的噪声FFT变换后第j个频点处的频域表达式；子段q内阵列接收信号的频域表达式为：

x_q(ω_j)＝A_q(ω_j,θ)s_q(ω_j)+n_q(ω_j)

其中x_q(ω_j)＝[X_{1_q}(ω_j),X_{2_q}(ω_j)，...，X_{M_q}(ω_j)]^T，

s_q(ω_j)＝[S_{1_q}(ω_j),S_{2_q}(ω_j)，...，S_{P_q}(ω_j)]^T，n_q(ω_j)＝[N_{1_q}(ω_j),N_{2_q}(ω_j)，...，N_{M_q}(ω_j)]^T，

θ＝[θ₁,θ₂,…,θ_P]^T，A_q(ω_j,θ)是M×P维关于波达角度和频率的矩阵，在第j个频点的表达式为：

3、计算信号子空间和噪声子空间：

第j个频点处阵列接收信号的协方差矩阵为：

4、定义对角变换矩阵，计算任意频点对应信号子空间：

为降低运算复杂度，在J个频点中等间距选取Z＝9个频点进行计算，并选取其中一个频点为参考频点，假定参考频点频率为ω₀，定义对角变换矩阵：

其中频率差Δω_z＝ω_z-ω₀，ω_z为选取的第z个频点频率。φ为待估计的波达角度。利用对角矩阵Φ(Δω_z,φ)将参考频点ω₀对应的信号子空间转化为任意频点对应的信号子空间，则第z个频点对应的信号子空间U_z(θ)为：

U_z(θ)＝Φ(Δω_z,φ)F₀z＝1,2，...,Z-1

F₀为参考频点ω₀对应的信号子空间。

5、为了便于说明，下面简单介绍一下原TOPS算法的处理方法：在原TOPS算法中，定义每一个频点处阵列导引矢量：

其正交投影矩阵为：

其中I为单位矩阵。投影后的信号子空间为：

其中U'_z(φ)为原TOPS算法第z个频点对应的投影后的信号子空间。

定义矩阵：

其中D'(φ)为原TOPS算法信号子空间与噪声子空间正交性测试矩阵，其中W_i为第i个频点对应的噪声子空间。当假定角度能匹配真实信号DOA，即φ∈θ时，矩阵D'(φ)不是一个行满秩矩阵，其秩会有所减少。可以通过矩阵D'(φ)的最小奇异值来判断该矩阵的缺秩程度，进而判断假定角度是不是入射信号的波达角度。因此，为了估计宽带信号DOA，需要先对D'(φ)进行奇异值分解，然后在阵列视角范围内进行空间角度搜索，波达角度的估计值位于下式谱峰位置：

其中σ_min(φ)为D(φ)的最小奇异值。

6、本发明中的改进TOPS算法在原有TOPS算法的基础上进行了两点改进：

6.1在原TOPS算法中，由于与信号子空间正交，得到的值会有所减小，即将信号子空间U_z投影到上会减弱信号子空间，通过可以看出，将噪声子空间W_z投影到上会加强噪声子空间，因而在一定程度上会导致虚假谱峰的出现。

为了解决这一问题，改进TOPS算法通过将信号子空间投影到信号增强区来增强信号子空间，减少虚假谱峰的出现。

投影后的信号子空间为：

其中为改进TOPS算法第z个频点对应的投影后的信号子空间。

定义矩阵：

其中D^Im(φ)为改进TOPS算法信号子空间与噪声子空间正交性测试矩阵，W_z为第z个频点对应的噪声子空间。并且当φ∈θ时，由于(I-P_az(φ))W_z＝0，则

可以看出改进TOPS算法第z个频点对应的投影后的信号子空间与噪声子空间W_z正交。

6.2在原TOPS算法中，利用D'(φ)的最小奇异值来衡量矩阵的缺秩程度，进而判断假定角度是不是入射信号的波达角度。即对D'(φ)进行奇异值分解，然后在阵列视角范围内进行空间角度搜索，波达角度的估计值位于下式谱峰位置：

然而奇异值分解算法的复杂程度较大，不利于工程实施。

为了解决这一问题，改进的TOPS算法利用D^Im(φ)与D^ImH(φ)点积向量的长度来衡量信号子空间与噪声子空间的正交性。在阵列视角范围内进行空间角度搜索，波达角度的估计值位于下式谱峰位置：

下面通过计算机仿真来验证本发明算法的性能。设ULA阵列的阵元数为8，进行仿真分析，信号中心频率为1GHz，带宽B＝200M，阵元间距为中心波长的一半，时域快拍为4096，均匀分为Q＝64段，分别对每段进行FFT计算，取7个频点进行分析仿真，

仿真参数如表1所示。

表1仿真参数表

图3(a)、图3(b)展示了不同信噪比下，TOPS算法和修正TOPS算法DOA估计的结果，图中三个信号源输入SNR相同。如图所示，虽然原来的TOPS算法在信号真实DOA处形成了谱峰，但是也形成了一些伪峰，而修正的mTOPS算法在真实的信号DOA处均形成了突出的谱峰，解决了原来TOPS算法虚假谱峰的问题。

图4为信源靠近时，TOPS算法与改进TOPS算法的性能对比图。如图4所示，两个宽带信号源的入射角分别以8°和9°，其输入SNR均为30dB，TOPS算法在估计时不仅形成了伪峰，且平均估计误差约0.3°左右，而修正TOPS算法的精度在此基础上有较大提升，能比较精确的估计两个信源相互靠近的宽带信号DOA信息。

为充分说明TOPS算法估计信号DOA的性能，采用均方根误差(RootMeanSquareError,RMSE)标准来量化分析DOA估计的精度，RMSE可定义为：

式中G为蒙特卡罗实验次数，表示第g次估计出来的目标空间角，θ_理论值表示理论目标空间角。图5为200次独立蒙特卡罗实验的平均值结果，展示了利用改进mTOPS算法与原来TOPS算法进行宽带信号DOA估计均方根误差与输入信噪比的关系，计算DOA时搜索精度为0.01°。由图可以看出，mTOPS算法较TOPS算法的DOA估计精度有较大程度的改善，并且在信噪比超过30dB时，其估计均方根误差接近于零。

图6(a)、图6(b)进一步分析了修正TOPS算法的计算复杂度。本次仿真中，以复数乘法的次数表示算法的计算复杂度，分别要对每个频点的协方差矩阵进行特征值分解，其计算量为O(M³)，其中M为阵元数。以Z表示计算的频点数，P为信源个数，TOPS算法需要(Z(M³+MP+PM²)+P(M-P)²Z²+(M-P)ZP²)次复乘运算，修正TOPS算法在进行搜索时不需要进行特征值分解，只需要进行矩阵相乘，因此其需要的复乘次数为(Z(M³+MP+PM²)+P(M-P)²Z²)。图6(a)、图6(b)分别分析了频点数为7时计算复杂度与阵元数之间关系以及阵元数为8时计算复杂度与频点数之间关系。由图可以看出，修正TOPS算法改善了TOPS算法的计算复杂度，尤其是当所用频点数较多时，修正TOPS算法能有效降低TOPS算法的计算复杂度。在实际系统中，当信源带宽较大时，为了更准确的估计每个频点的信号子空间，用于计算的频点数较多，为了满足算法所需的自由度，阵元数一定要超过信源数，此时，修正TOPS算法能有效降低算法的计算量，提高系统的运算效率。

本技术领域技术人员可以理解的是，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims

1.一种宽带阵列波达方向估计方法，其特征在于，包括如下步骤：

<mrow> <msub> <mi>X</mi> <mrow> <mi>m</mi> <mo>_</mo> <mi>q</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>&omega;</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>P</mi> </munderover> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>j&omega;</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mfrac> <mi>d</mi> <mi>c</mi> </mfrac> <msub> <mi>sin&theta;</mi> <mi>p</mi> </msub> </mrow> </msup> <msub> <mi>S</mi> <mrow> <mi>p</mi> <mo>_</mo> <mi>q</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>&omega;</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>N</mi> <mrow> <mi>m</mi> <mo>_</mo> <mi>q</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>&omega;</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> </mrow>

其中，ω_j为阵元接收信号FFT变换后第j个频点处的信号频率，ω_j∈[ω_L,ω_H]，j＝1,2...,J为子段q内FFT变换的频点数，ω_L、ω_H分别为子段q内第m个阵元接收信号的最低频率与最高频率，且ω_L＝ω₁，ω_H＝ω_J；S_{p_q}(ω_j)为子段q内第p个信源在第一个阵元处的接收信号FFT变换后第j个频点处的频域表达式，N_{m_p}(ω_j)为子段q内第m个阵元接收的噪声FFT变换后第j个频点处的频域表达式；

其中，x_q(ω_j)＝[X_{1_q}(ω_j),X_{2_q}(ω_j),...,X_{M_q}(ω_j)]^T，s_q(ω_j)＝[S_{1_q}(ω_j),S_{2_q}(ω_j),...,S_{P_q}(ω_j)]^T，n_q(ω_j)＝[N_{1_q}(ω_j),N_{2_q}(ω_j),...,N_{M_q}(ω_j)]^T，θ＝[θ₁,θ₂,…,θ_P]^T，A_q(ω_j,θ)是M×P维关于波达角度和频率的矩阵，在第j个频点的表达式为

步骤3)，计算信号子空间和噪声子空间：

第j个频点处阵列接收信号的协方差矩阵为j＝1,2,...,J；

步骤4.1)，在J个频点中等间距选取Z＝9个频点进行计算，并选取其中一个频点为参考频点，令参考频点频率为ω₀，定义对角变换矩阵z＝1,2，...,Z-1，ω_z为选取的第z个频点频率，频率差Δω_z＝ω_z-ω₀，φ为待估计的波达角度；

步骤4.2)，利用对角矩阵Φ(Δω_z,φ)将参考频点ω₀对应的信号子空间转化为任意频点对应的信号子空间，则第z个频点对应的信号子空间U_z(θ)为：

U_z(θ)＝Φ(Δω_z,φ)F₀

式中，z＝1,2，...,Z-1，F₀为参考频点ω₀对应的信号子空间；

步骤5)，计算信号子空间与噪声子空间正交性测试矩阵：

步骤5.1)，定义第z个频点处的阵列导引矢量z＝1,2，...,Z-1，其正交投影矩阵I为单位矩阵；

2.根据权利要求1所述的宽带阵列波达方向估计方法，其特征在于，所述步骤2)中，