基于延迟快拍的均匀线列阵高分辨波达方向估计方法与流程

本发明涉及阵列信号处理领域，特别涉及一种基于延迟快拍的均匀线列阵高分辨波达方向估计方法。

背景技术：

阵列信号处理广泛应用于移动通信、雷达、声呐、导航、地震探测和医学成像等领域。波达方向估计是阵列信号处理的关键技术之一。具体来说，波达方向估计是指，对阵列传感器接收信号进行分析与处理从而估计信号到达方向的过程。

已有的波达方向估计方法主要分为两类：一是以波束形成为代表的常规波达方向估计方法；二是以多重信号分类(multiplesignalclassification,music)为代表的高分辨波达方向估计方法。其中，波束形成方法的优点是计算量小、实现简单，缺点是当空间中存在两个或两个以上位置相邻的信源时，波束形成方法不能对它们进行有效的分辨。多重信号分类方法的方向分辨能力强于波束形成方法。但是，当空间中信源位置进一步靠近时，不论是波束形成方法还是多重信号分类方法都不能准确地估计出信号波达方向，而是将多个位置相邻的信源看作一个信源。因此，研究一种方向分辨能力更高的波达方向估计方法具有重要的实用价值。

技术实现要素：

为了有效解决对于邻近信源，已有波达方向估计方法性能显著恶化的问题，本发明提供了一种基于延迟快拍的均匀线列阵高分辨波达方向估计方法。在同等信噪比和同等快拍数的情况下，本发明方法的波达方向估计精度明显优于已有波达方向估计方法。

本发明采用以下技术方案实现上述目的。基于延迟快拍的均匀线列阵高分辨波达方向估计方法，其特征在于，其步骤如下：

1)获取并存储阵列快拍：

首先，对所有阵元的接收信号进行采样，获取t0时刻的阵列快拍数据，即：

{x(1,t0),x(2,t0),...,x(n,t0)}；

其中：括号中的第一个数字表示阵元序号，n表示所有阵元的个数；括号中的第二个符号表示采样时刻；x(n,t0)表示第n个阵元在t0时刻的信号采样值；

然后，将t0时刻的阵列快拍数据保存在存储器中；

其次，获取并存储在t0+ts时刻的阵列快拍数据，即：

{x(1,t0+ts),x(2,t0+ts),...,x(n,t0+ts)}；

其中：ts表示采样周期，其大小等于采样频率fs的倒数，即ts＝1/fs；采用类似的方式，获取并存储在t0+its时刻的阵列快拍数据，其中i＝2,3,…k,其中，k表示最大的阵列快拍次数；

2)构造延迟快拍：

首先，利用存储器中保存的各个时刻的阵列快拍数据，构造t0时刻的延迟快拍数据，即：

{x(1,t0),x(2,t0+αts),...,x(n,t0+α(n-1)ts)}；

其中，α表示延迟快拍因子，延迟快拍因子应取整数，其具体大小根据实际应用情况确定；

将t0时刻的延迟快拍数据作为数据矩阵a的第一列，即：

a(:,1)＝[x(1,t0),x(2,t0+αts),...,x(n,t0+α(n-1)ts)]^t；

其中，符号a(:,1)表示数据矩阵a的第一列，符号[·]^t表示转置操作；

然后，构造t0+ts时刻的延迟快拍数据，即：

{x(1,t0+ts),x(2,t0+αts+ts),...,x(n,t0+α(n-1)ts+ts)}；

并将t0+ts时刻的延迟快拍数据作为数据矩阵a的第二列，即：

a(:,2)＝[x(1,t0+ts),x(2,t0+αts+ts),...,x(n,t0+α(n-1)ts+ts)]^t；

其中，符号a(:,2)表示数据矩阵a的第二列；以此类推，将t0+its时刻的延迟快拍数据作为

数据矩阵a的第i列，即：

a(:,i)＝[x(1,t0+its),x(2,t0+αts+its),...,x(n,t0+α(n-1)ts+its)]^t；

其中，符号a(:,i)表示数据矩阵a的第i列；这里，i＝2,3,…,m,其中，m既是最大的延迟快拍个数，也是数据矩阵a的列的个数，其具体大小由实际需要的延迟快拍数决定；

3)计算信号协方差矩阵：

根据数据矩阵a计算信号协方差矩阵r，即：

4)对信号协方差矩阵进行特征分解：

对信号协方差矩阵r进行特征分解，即：

r＝qλq^t；

其中，矩阵q＝[q1,q2,...,qn]，列向量qi表示信号协方差矩阵r的第i个特征向量；矩阵λ＝diag(λ1,λ2,...,λn)，符号diag(·)表示矩阵λ为对角矩阵，矩阵λ的对角元素分别为λ1,λ2,...,λn，数值λi表示信号协方差矩阵r的第i个特征值；矩阵λ的对角元素按照从大到小排列，即λ1≥λ2≥...≥λn；

5)构造噪声子空间：

由信号协方差矩阵r的特征向量ql+1,ql+2,...,qn，构造噪声子空间rn，其中，l表示信源个数，其具体大小根据实际应用情况确定；噪声子空间rn的具体构造如下：

rn＝[ql+1,ql+2,...,qn]；

6)计算空间功率谱：

根据噪声子空间rn和阵列导向向量a(θ)计算空间功率谱p(θ)，即：

其中，符号[·]^h表示共轭转置操作；均匀线列阵对应的导向向量如下：

a(θ)＝[1,e^{j(2π/λ)dsinθ},...,e^{j(2π/λ)(n-1)dsinθ}]^t；

其中，符号d表示阵元间距，单位为：米；符号λ表示信号的波长，单位为：米；符号θ表示信号的波达方向，单位为：度；

7)搜索谱峰，确定映射波达方向：

在波达方向范围[-90°,90°]内，观察空间功率谱p(θ)的大小，搜索空间功率谱p(θ)的l个谱峰，将所有l个谱峰对应的角度作为信号的映射波达方向，并记为

8)计算信号的估计波达方向：

根据信号的映射波达方向计算信号的实际波达方向，即：

其中，θ1,θ2,...,θl为信号的l个实际波达方向，arcsin(·)表示反正弦函数，c表示信号在介质中的传播速度，单位为：米/秒。

进一步，所述存储器的存储空间大于或等于n²×d字节，其中：d表示存储单个采样数据所需的字节数，d的具体大小根据实际应用情况确定。

进一步，所述采样频率的大小由实际目标源信号的频率决定。

本发明的有益效果是：在同等信噪比及同等快拍数情况下，相比于已有的波达方向估计方法，本发明方法的空间分辨能力更高，能够取得更高的波达方向估计精度。另外，本发明方法的计算复杂度与多重信号分类方法相当，因此本发明方法在实施起来比较方便。

附图说明

图1为本发明基于延迟快拍的均匀线列阵高分辨波达方向估计方法的流程图；

图2为延迟快拍因子α＝1时均匀线列阵延迟快拍的示意图；

图3为仿真实验中本发明方法和music方法的空间功率谱；

图4为图3中music方法空间功率谱的局部放大图；

图5为图3中本发明方法空间功率谱的局部放大图；

图6为仿真实验中实际波达方向与映射波达方向之间的对应曲线。

具体实施方式

下面结合附图对本发明进行详细说明。基于延迟快拍的均匀线列阵高分辨波达方向估计方法，首先，对所有阵元的接收信号进行采样，将获得的阵列快拍信号存储到存储器中；其次，按照一定的顺序选取存储器中的阵列快拍信号，构造延迟快拍，并由延迟快拍得到数据矩阵；接着，由数据矩阵计算得到信号协方差矩阵，对信号协方差矩阵进行特征分解得到特征值和特征向量，并由较小特征值对应的特征向量构成噪声子空间；然后，根据噪声子空间和阵列导向向量计算空间功率谱，并搜索空间功率谱的谱峰，确定映射波达方向；最后，按照一定的对应关系将映射波达方向转换为实际波达方向。

图1是本发明所提供的基于延迟快拍的均匀线列阵高分辨波达方向估计方法的流程图。该流程图的第一步是获取并存储阵列快拍101，首先，对所有阵元的接收信号进行采样，然后，将获得的阵列快拍信号存储到存储器中。第二步是构造延迟快拍102，即按照一定的顺序选取存储器中的阵列快拍信号形成延迟快拍，并由延迟快拍计算得到数据矩阵a。第三步是计算信号协方差矩阵103，即由数据矩阵a计算信号协方差矩阵r。第四步是信号协方差矩阵的特征分解104，即对信号协方差矩阵r进行特征分解，得到特征向量q1,q2,...,qn及特征值λ1,λ2,...,λn。第五步是构造噪声子空间105，即由较小特征值对应的特征向量ql+1,ql+2,...,qn构成噪声子空间rn。第六步是计算空间功率谱106，即由噪声子空间rn和阵列导向向量a(θ)计算空间功率谱p(θ)。第七步是搜索谱峰，确定映射波达方向107，即搜索空间功率谱p(θ)的谱峰，并将谱峰所对应的角度作为映射波达方向。第八步，计算信号的估计波达方向108，即按照一定的对应关系将映射波达方向转换为实际波达方向。

图2是延迟快拍因子α＝1时均匀线列阵延迟快拍的示意图。其中，第一行的圆点表示阵元1、阵元2、…、以及阵元n在起始时刻t0时的阵列快拍数据，在图2中记为“快拍1”。第二行的圆点表示在时刻t0+ts时的阵列快拍数据，在图2中记为“快拍2”。以此类推，第k行的圆点表示在时刻t0+(k-1)ts时的阵列快拍数据，记为“快拍k”。在图2中从上往下看，第一条虚斜线经过的圆点构成t0时刻的延迟快拍数据，第二条虚斜线经过的圆点构成t0+ts时刻的延迟快拍数据，以此类推，第i条虚斜线经过的圆点构成t0+(i-1)ts时刻的延迟快拍数据。

图3至图6是利用本发明方法在具体实施方式下进行仿真实验的处理结果，其中仿真实验的基本参数设置如下：

均匀线列阵的阵元数n＝13，阵元间距d＝0.025米，阵列接收信号由正弦信号和高斯白噪声组成，其中正弦信号的个数l＝2。两个正弦信号的频率分别为f1＝30000hz和f1＝30100hz，两个正弦信号的真实波达方向分别为度和度，两个正弦信号的信噪比分别为10db和20db。信号采样频率为fs＝260000hz，延迟快拍数为k＝20000，延迟快拍因子α＝4，信号在介质中的传播速度c＝1500米/秒。

图3是仿真实验中本发明方法和music方法的空间功率谱。其中，实线表示music方法的空间功率谱，虚线表示本发明方法的空间功率谱。图4是图3中music方法空间功率谱的局部放大图。图5是图3中本发明方法空间功率谱的局部放大图。由图3及图4可以发现：(1)music方法的谱峰个数与信源个数不匹配，music方法的空间功率谱仅存在一个明显的谱峰，而实际的信源个数为2个；(2)music方法的波达方法估计精度较低，music方法空间功率谱的谱峰位于1.66度，与真实波达方向(度和度)相比均存在一定的误差。由图3及图5可以发现：(1)本发明方法的谱峰个数与信源个数相一致，本发明方法的空间功率谱存在两个明显的谱峰，与实际的信源个数一致；(2)本发明方法的波达方法估计精度更高，本发明方法空间功率谱的谱峰分别位于69.87度和72.65度，即映射波达方向度,度。将映射波达方向进行转换，得到本发明方法的波达方向估计结果，即θ1＝0.907446度和θ2＝1.800759度，波达方向估计结果非常接近真实波达方向(度和度)。

图6为仿真实验中实际波达方向与映射波达方向之间的对应曲线。在图6中，横轴表示实际波达方向，纵轴表示映射波达方向。图6中的曲线是当延迟快拍因子α＝4，信号采样频率为fs＝260000hz，信号在介质中的传播速度c＝1500米/秒，阵元间距d＝0.025米时，实际波达方向与映射波达方向之间的对应关系。当映射波达方向分别为度,度时，参照曲线可得实际波达方向分别为θ1＝0.907446度和θ2＝1.800759度。