CN109633522A - 基于改进的music算法的波达方向估计方法 - Google Patents

Abstract

本发明属于通信技术领域，公开了基于改进的MUSIC算法的波达方向估计方法。该方法包括：建立对称阵列模型，对各天线阵元接收的射频信号进行采样，正交下变频和滤波得到复信号，对复信号构成的阵列输出矩阵做线性变换得到复矩阵，进而处理得到实对称矩阵，对实对称矩阵进行特征值分解，构建噪声子空间矩阵，构造阵列方向矩阵，利用噪声子空间和阵列方向矩阵构造空间谱函数，进而进行谱峰局部极值搜索，得到函数极大值对应的角度值即为波达方向。本发明能够在不损失算法性能的情况下将特征值分解由复矩阵运算转换为实矩阵运算，并利用实方向矩阵进行峰值搜索，从而能有效地估计出信号的波达方向，提高了波达方向估计的准确性并且降低了计算的复杂度。

Description

基于改进的MUSIC算法的波达方向估计方法

技术领域

本发明涉及通信技术领域，尤其涉及基于改进的MUSIC算法的波达方向估计方法。

背景技术

在无线电通信、导航、雷达、声纳、电子侦察等领域，获得信号的来波方向极其重要。在军事无线通信领域，获得高分辨率的来波方向不仅有助于对我方目标的导航、跟踪，还可以有助于我方对敌方目标进行准确定位以及打击。在民用无线通信领域，准确的来波方向估计能够提升系统容量和减少多址干扰等。

空间谱表示信号在各个方向上的能量分布，若能够估计出空间谱就可以得到信号源的波达方向，所以空间谱估计也常称为波达方向(英文全称：Direction Of Arrival，英文缩写：DOA)估计，DOA估计技术的关键在于利用空间不同位置的天线阵列，接收来自不同方向的多个信号源发出的信号，运用现代信号处理技术估计出信号源的来波方向。在多种空间谱估计算法中，多重信号分类(英文全称：Multiple Signal Classification，英文缩写：MUSIC)算法是一种基于子空间分解的算法，其特点为测向误差小、精度高、分辨率高、灵敏度高，可以对入射信号数目、到达方向及波形的强度做出无偏估计。因此研究MUSIC算法对于促进更广泛的应用有着十分重要的意义。

但是，传统的MUSIC算法在协方差矩阵计算及特征分解时涉及大量的复数运算，导致算法的复杂度高，不利于硬件实现。因此，如何在降低计算复杂度的同时又不损失算法性能就显得尤为重要。

发明内容

本发明的实施例提供基于改进的MUSIC算法的波达方向估计方法，能够在不损失算法性能的情况下将特征值分解由复矩阵运算转换为实矩阵运算，并利用实方向矩阵进行峰值搜索，从而能有效地估计出信号的波达方向。此发明提高了波达方向估计的准确性并且降低了计算的复杂度。

为达到上述目的，本发明的实施例采用如下技术方案：

步骤1、建立对称阵列模型，天线阵列为均匀线阵或者平面均匀圆阵，阵元个数为偶数；当天线阵元为均匀线阵时，以阵列的对称轴为基准，向右将阵列对称轴右侧天线阵元依次编号为向左将阵列对称轴左侧天线阵元依次编号为当天线阵元为平面均匀圆阵时，以其中任一阵元为起点，按顺时针方向依次编号为1，2...M，阵元间隔为d；

步骤2、对各天线阵元接收的射频模拟信号进行K次采样，其中，第m个阵元天线的第k次采样得到的采样信号为s_m(k)，m＝1，2，…M，k＝1，2，…K；K为快拍数，K＞2f_s/Δf，f_s为采样频率，Δf为需要分辨的最小频率，M为天线阵元个数，M为偶数；

步骤3、对各天线阵元的K个采样信号做I、Q正交下变频和滤波，得到各天线阵元的K个复信号，并对第一阵元的K个复信号做快速傅里叶变换，判断快速傅里叶变换的最大值是否超过预设的门限值，如果是，则执行步骤4；

步骤4、利用各天线阵元的K个复信号构成阵列输出矩阵X，对阵列输出矩阵X做线性变换得到复矩阵Y，利用复矩阵Y进一步计算得到复共轭对称矩阵R，利用复共轭对称矩阵R进一步计算得到实对称矩阵R_X；

步骤5、对实对称矩阵R_X进行特征值分解，得到实对称矩阵R_X的特征值和对应的特征向量，并将特征值按照从大到小的顺序进行排序；

步骤6、基于最小描述长度MDL准则，利用排序后的实对称矩阵R_X的特征值估计信号源个数，进而根据信号源个数，利用实对称矩阵R_X的特征向量得到噪声子空间矩阵；

步骤7、基于对称阵列模型构造阵列方向矩阵，利用噪声子空间和阵列方向矩阵构造空间谱函数，进而利用空间谱函数进行谱峰局部极值搜索，得到空间谱函数极大值对应的角度值，角度值即为信号源的入射方向。

本发明的实施例提供基于改进的MUSIC算法的波达方向估计方法，首先建立天线的对称阵列模型，通过对中心对称阵列结构的数学分析，利用采样数据及其共轭信息得到了一种线性变换方法，能够在不损失算法性能的情况下将特征值分解由复矩阵运算转换为实矩阵运算，并利用实方向矩阵进行峰值搜索，从而能有效地估计出信号的波达方向。此发明提高了波达方向估计的准确性并且降低了计算的复杂度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的基于改进的MUSIC算法的波达方向估计方法的流程示意图；

图2为四阵元均匀直线阵的阵列结构示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1为本发明实施例提供的基于改进的MUSIC算法的波达方向估计方法的流程示意图，如图1所示，

步骤1、建立对称阵列模型，天线阵列为均匀线阵或者平面均匀圆阵，阵元个数为偶数；当天线阵元为均匀线阵时，以阵列的对称轴为基准，向右将阵列对称轴右侧天线阵元依次编号为向左将阵列对称轴左侧天线阵元依次编号为当天线阵元为平面均匀圆阵时，以其中任一阵元为起点，按顺时针方向依次编号为1，2...M，阵元间隔为d；

步骤2、对各天线阵元接收的射频模拟信号进行K次采样，其中，第m个阵元天线的第k次采样得到的采样信号为s_m(k)，m＝1，2，…M，k＝1，2，…K；K为快拍数，K＞2f_s/Δf，f_s为采样频率，Δf为需要分辨的最小频率，M为天线阵元个数，M为偶数；

步骤3、对各天线阵元的K个采样信号做I、Q正交下变频和滤波，得到各天线阵元的K个复信号，并对第一阵元的K个复信号做快速傅里叶变换，判断快速傅里叶变换的最大值是否超过预设的门限值，如果是，则执行步骤4；

进一步的，利用公式对第一阵元的K个复信号做快速傅里叶变换；

其中，第m个天线阵元第k次采样对应的复信号为x_m(k)＝I_m(k)+jQ_m(k)，m＝1，2，…M，k＝1，2，…K；N为快速傅里叶变换的点数，

步骤4、利用各天线阵元的K个复信号构成阵列输出矩阵X，对阵列输出矩阵X做线性变换得到复矩阵Y，利用复矩阵Y进一步计算得到复共轭对称矩阵R，利用复共轭对称矩阵R进一步计算得到实对称矩阵R_X；

进一步的，步骤4具体包括：

(4.1)利用各天线阵元的K个复信号构成阵列输出矩阵X，X＝[X₁(k)，X₂(k)，...，X_m(k)，...，X_M(k)]^T，X_m(k)＝[x_m(1)，x_m(2)，...x_m(k)...，x_m(K)]，m＝1，2，…M；

(4.2)构造线性变换矩阵根据线性变换矩阵T和阵列输出矩阵X，计算得到复矩阵Y＝TX；

其中，M为表示阵元个数，M为偶数；

(4.3)利用复矩阵Y，计算得到复共轭对称矩阵R＝Y×Y^H，其中，上标H表示矩阵的共轭转置；

(4.4)取复共轭对称矩阵R的实部，得到实对称矩阵R_X＝Re{R}，其中，Re{·}表示取矩阵的实部。

步骤5、对实对称矩阵R_X进行特征值分解，得到实对称矩阵R_X的特征值和对应的特征向量，并将特征值按照从大到小的顺序进行排序；

进一步的，步骤5具体包括：

(5.1)利用预设公式：R_X＝UAU^H，对实对称矩阵R_X进行特征值分解；

其中，Λ为对角矩阵，其对角线上的元素为实对称矩阵R_X的特征值，U是R_X的特征矢量矩阵；上标H表示矩阵的共轭转置；

(5.2)将R_X的特征值按大小顺序排列，得到λ₁≥λ₂≥…≥λ_M。

步骤6、基于最小描述长度MDL准则，利用排序后的实对称矩阵R_X的特征值估计信号源个数，进而根据信号源个数，利用实对称矩阵R_X的特征向量得到噪声子空间矩阵。

进一步的，步骤6具体包括：

(6.1)基于最小描述长度MDL准则，利用特征值构造信号源数目估计模型：

其中，D的取值是0，1，…，M-1中的一个。

(6.2)利用信号源数目估计模型，计算得到空间信号源数目的估计值

其中，arg{min MDL(D)}表示MDL(D)取最小值时，D的取值；

(6.3)利用实对称矩阵R_X的个较小特征值对应的特征向量构成噪声子空间矩阵噪声子空间矩阵U_n的维度为

步骤7、基于对称阵列模型构造阵列方向矩阵，利用噪声子空间和阵列方向矩阵构造空间谱函数，进而利用空间谱函数进行谱峰局部极值搜索，得到空间谱函数极大值对应的角度值，角度值即为信号源的入射方向。

优选的，步骤7具体包括：

(7.1)基于对称阵列模型构造阵列方向矩阵第n个列矢量为，进一步计算得到阵列的实方向矩阵

(7.2)利用实方向矩阵的各个列矢量和噪声子空间矩阵U_n构造空间谱函数：

其中，θ为信号的入射角度；

(7.3)计算P_MUSIC(θ)的个极值点，进一步得到个极值点对应的角度值角度值即为信号源的入射方向。

参照图2，以四阵元均匀直线阵(M＝4)为例，对线性变换的计算过程进行原理性说明：

为了符合阵列在结构上的对称性，将阵元1和阵元2组成子阵1，阵元3和阵元4组成子阵2，信号入射角度用θ表示，则子阵1的方向矩阵A₁和子阵2的方向矩阵A₂可以分别表示为

其中，阵元间距λ表示入射信号波长。

进一步得到子阵1的输出矩阵X₁和子阵2的输出矩阵X₂：

其中，S(k)为时域连续K个采样快拍的入射信号构成的矩阵，N₁(k)为子阵1的噪声矩阵，N₂(k)和子阵2噪声矩阵。

由于两个子阵的方向矩阵A₁和A₂是互为共轭对称的。于是构造一个线性变换矩阵T：

对X(k)作线性变换：

则有，为一个实矩阵，其中：

N(k)为阵列的噪声矩阵。

本发明的实施例提供基于改进的MUSIC算法的波达方向估计方法，首先建立天线的对称阵列模型，通过对中心对称阵列结构的数学分析，利用采样数据及其共轭信息得到了一种线性变换方法，能够在不损失算法性能的情况下将特征值分解由复矩阵运算转换为实矩阵运算，并利用实方向矩阵进行峰值搜索，从而能有效地估计出信号的波达方向。此发明提高了波达方向估计的准确性并且降低了计算的复杂度。

本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。

以上，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。

Claims (6)

1.一种基于改进的MUSIC算法的波达方向估计方法，其特征在于，包括以下步骤：

步骤1、建立对称阵列模型，天线阵列为均匀线阵或者平面均匀圆阵，阵元个数为偶数；当天线阵元为均匀线阵时，以阵列的对称轴为基准，向右将阵列对称轴右侧天线阵元依次编号为向左将阵列对称轴左侧天线阵元依次编号为当天线阵元为平面均匀圆阵时，以其中任一阵元为起点，按顺时针方向依次编号为1，2...M，阵元间隔为d；

步骤2、对各天线阵元接收的射频模拟信号进行K次采样，其中，第m个阵元天线的第k次采样得到的采样信号为s_m(k)，m＝1，2，…M，k＝1，2，…K；K为快拍数，K＞2f_s/Δf，f_s为采样频率，Δf为需要分辨的最小频率，M为天线阵元个数，M为偶数；

步骤3、对所述各天线阵元的K个采样信号做I、Q正交下变频和滤波，得到所述各天线阵元的K个复信号，并对第一阵元的K个复信号做快速傅里叶变换，判断所述快速傅里叶变换的最大值是否超过预设的门限值，如果是，则执行步骤4；

步骤4、利用所述各天线阵元的K个复信号构成阵列输出矩阵X，对阵列输出矩阵X做线性变换得到复矩阵Y，利用复矩阵Y进一步计算得到复共轭对称矩阵R，利用复共轭对称矩阵R进一步计算得到实对称矩阵R_X；

步骤5、对实对称矩阵R_X进行特征值分解，得到实对称矩阵R_X的特征值和对应的特征向量，并将所述特征值按照从大到小的顺序进行排序；

步骤6、基于最小描述长度MDL准则，利用排序后的实对称矩阵R_X的特征值估计信号源个数，进而根据信号源个数，利用实对称矩阵R_X的特征向量得到噪声子空间矩阵；

步骤7、基于对称阵列模型构造阵列方向矩阵，利用噪声子空间和阵列方向矩阵构造空间谱函数，进而利用空间谱函数进行谱峰局部极值搜索，得到所述空间谱函数极大值对应的角度值，所述角度值即为所述信号源的入射方向。

2.根据权利要求1所述的方法，其特征在于，在步骤3中，利用公式对第一阵元的K个复信号做快速傅里叶变换；

其中，第m个天线阵元第k次采样对应的复信号为x_m(k)＝I_m(k)+jQ_m(k)，m＝1，2，…M，k＝1，2，…K；N为快速傅里叶变换的点数，

3.根据权利要求2所述的方法，其特征在于，所述步骤4包括：

(4.1)利用所述各天线阵元的K个复信号构成阵列输出矩阵X，X＝[X₁(k)，X₂(k),...，X_m(k),...，X_M(k)]^T，X_m(k)＝[x_m(1)，x_m(2),...x_m(k)...，x_m(K)]，m＝1，2，…M；

(4.2)构造线性变换矩阵根据线性变换矩阵T和阵列输出矩阵X，计算得到复矩阵Y＝TX；

其中，M为表示阵元个数，M为偶数；

(4.3)利用复矩阵Y，计算得到复共轭对称矩阵R＝Y×Y^H，其中，上标H表示矩阵的共轭转置；

(4.4)取所述复共轭对称矩阵R的实部，得到实对称矩阵R_X＝Re{R}，其中，Re{·}表示取矩阵的实部。

4.根据权利要求3所述的方法，其特征在于，所述步骤5中具体包括：

(5.1)利用预设公式：R_X＝UAU^H，对实对称矩阵R_X进行特征值分解；

其中，Λ为对角矩阵，其对角线上的元素为实对称矩阵R_X的特征值，U是R_X的特征矢量矩阵；上标H表示矩阵的共轭转置；

(5.2)将R_X的特征值按大小顺序排列，得到λ₁≥λ₂≥…≥λ_M。

5.根据权利要求4所述的方法，其特征在于，所述步骤6具体为：

(6.1)基于最小描述长度MDL准则，利用特征值构造信号源数目估计模型：

其中，D的取值是0，1，…，M-1中的一个；

(6.2)利用信号源数目估计模型，计算得到空间信号源数目的估计值

其中，arg{minMDL(D)}表示MDL(D)取最小值时，D的取值；

(6.3)利用实对称矩阵R_X的个较小特征值对应的特征向量构成噪声子空间矩阵噪声子空间矩阵U_n的维度为

6.根据权利要求5所述的方法，其特征在于，所述步骤7具体包括：

(7.1)基于对称阵列模型构造阵列方向矩阵第n个列矢量为，进一步计算得到阵列的实方向矩阵

(7.2)利用实方向矩阵的各个列矢量和噪声子空间矩阵U_n构造空间谱函数：

其中，θ为信号的入射角度；

(7.3)计算P_MUSIC(θ)的个极值点，进一步得到个极值点对应的角度值所述角度值即为所述信号源的入射方向。