CN106970349A - 一种基于改进的music算法的ads‑b信号波达方向估计方法 - Google Patents

Abstract

本发明公开了一种基于改进的MUSIC算法的ADS‑B信号波达方向估计方法，本发明包括以下几个步骤：一、建立线型阵列的ADS‑B信号采集模型，二、估计ADS‑B信号的空间谱函数，三、基于求二阶导数的方法处理ADS‑B信号的空间谱函数及对二阶导数空间谱函数进行谱峰搜索。本发明用于民航飞机ADS‑B信号人为干扰抑制及运行安全保障。

Description

一种基于改进的MUSIC算法的ADS-B信号波达方向估计方法

技术领域

本发明属于民用航空器ADS-B监视及运行安全保障领域，特别是涉及一种基于改进的MUSIC算法的ADS-B信号波达方向估计方法与可行性验证方法。

背景技术

为了应对日益增长的航空需求和空中交通管制相对落后的矛盾，国际民航组织提出一种新航行系统下的广播式自动相关监视技术ADS-B。ADS-B是一种基于全球卫星定位系统和空-空、地-空数据链通信的航空器运行监视技术。ADS-B和传统雷达监视系统相比可以获得更高精度的目标位置和高度信息，突破其无法覆盖荒漠远洋等地区、数据精度受限的问题。通过为飞行员提供更加快捷的飞行数据传输服务，使用了ADS-B可以优化航路设置，减小间隔标准，使空域得到更充分的利用。但是ADS-B信号因为其非加密的编码方式，在人为干扰面前非常脆弱，极易受到多种人为干扰的攻击。如压制性干扰、欺骗性干扰和多径干扰等。

在干扰抑制方面常用的方法种类有很多，最常见的是一种基于阵列天线的波束形成方法，调整天线阵列权重，使阵列天线对准期望信号的来波方向，从而达到抑制人为干扰的目的。波束形成中最重要的一个环节是估计ADS-B信号的来波方向。经典的估计来波方向的方法有最小方差法Capon算法，多重信号分类(Multi Signal Classification，简称MUSIC)算法和旋转不变子空间信号参数估计(Estimation of Signal Parameters viaRotational Invariance Techniques，简称ESPRIT)算法等。考虑到Capon算法需要对矩阵求逆运算，大大增加了算法的计算量，而MUSIC和ESPRIT方法使谱估计突破瑞利限约束，同时综合考虑本例中采用现行阵列天线的实际情况，本例中采用MUSIC算法。

MUSIC算法是基于空间信号中信号子空间和噪声子空间的正交特性，利用阵列协方差矩阵的特征值分解，对信号的波达方向进行超分辨估计的一种方法。但是MUSIC算法在算法复杂度和估计精度上也存在局限，针对此，近些年也相继出现了很多基于MUSIC算法的改进型算法。如用多项式求根来代替MUSIC算法中的谱搜索的求根-MUSIC算法，在低信噪比的情况下，提升了算法的分辨率；利用接收信号的谱相干性和空间相干性的信号选择定向算法的循环MUSIC算法，这种算法不受信号数小于阵元数这个条件的约束；还有基于空间平滑技术的空间平滑MUSIC算法，旨在解决空间信号相关造成的算法估计性能下降的问题。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供一种改进型的MUSIC算法的ADS-B信号波达方向估计方法，用于解决现有技术中MUSIC算法分辨率低和性能不佳的技术问题。

技术方案：为实现上述目的，本发明采用的技术方案为：步骤一、建立线型阵列的ADS-B信号采集模型

假设接收信号满足窄带条件，即信号经过阵列长度所需要的时间应远远小于信号的相干时间，信号包络在天线阵列的传播时间内变化不大。为简化，假定信源和天线阵列是在同一平面内，并且入射到天线阵列的入射波为平面波。

以来波方向入射M根天线、阵元间距为d的均匀线阵，入射信号波长为λ，入射信号的信源总数为K个，第k个信源的到达角为θ_k，则天线阵列的响应矢量a(θ_k)为

定义方向矩阵为A

A＝[a(θ₁)，a(θ₂)，...a(θ_k),...a(θ_K)] (2)

即

用X(t)来表示阵元接收到的信号，则远场窄带信号线型阵列采集模型可以表示为

X(t)＝AS(t)+N(t) (4)

其中，S(t)为信号向量，N(t)为噪声向量，且有

S(t)＝[s₁(t)，s₂(t)，...s_k(t)，...s_K(t)]^T (5)

N(t)＝[n₁(t)，n₂(t)，...n_m(t)，...,n_M(t)]^T (6)

其中，s_k(t)为阵列接收到的第k个信源信号，n_m(t)为第m个阵元接收到的加性白噪声；

步骤二、估计ADS-B信号的空间谱函数；

本发明是基于一种改进的MUSIC算法对ADS-B信号的来波方向进行估计，具体步骤是在传统MUSIC算法估计出的空间谱函数的基础之上，对空间谱函数进行处理，所以首先需要得到ADS-B信号的空间谱函数。基于传统MUSIC算法的理论，利用信号子空间和噪声子空间的正交性，先求出接收信号的协方差矩阵，对协方差矩阵进行特征值分解，再对特征值进行排序，分辨出噪声子空间和信号子空间，最后得出空间谱函数。

具体实现分为以下几个步骤

步骤A、根据矩阵特征分解的理论，对阵列协方差矩阵进行特征分解。

首先，无噪声情况下的协方差矩阵R_s表示为

R_s＝AR_sA^H (7)

那么有噪声存在的情况下的协方差矩阵R_X表示为

R_X＝AR_sA^H+σ²I (8)

其中：

上标H表示共轭转置；I表示单位矩阵；σ²表示噪声功率，由于σ²>0，R_X为满秩矩阵，所以R_X有M个正特征值λ₁,λ₂,...,λ_M，分别对应于M个特征向量v₁,v₂,...,v_M，又由于R_X是Hermite矩阵，所以R_X的各特征向量是相互正交的。

由于入射信号的信源数为K个，则R_X的特征值中与信号有关的特征值只有K个，且矩阵AR_sA^H的各特征值分别与σ²求和所得值分别等于这K个与信号有关的特征值，R_X的特征值中其余的M-K个特征值为σ²，也就是说，σ²是R_X的最小特征值。

同理，特征向量v₁,v₂,...,v_M中，也只有K个是与信号有关的，另外M-K个是与噪声有关的。

步骤B、对阵列协方差矩阵进行特征分解，得到相互正交的信号子空间和噪声子空间：

将矩阵R_X的特征值进行从小到大的排序，即

λ₁≥λ₂≥...≥λ_M>0 (9)

其中K个较大的特征值对应于信号，M-K个较小的特征值对应于噪声。同理，矩阵R_X的属于这些特征值的特征向量也分别对应于信号和噪声，因此，可以把R_X的特征值(特征向量)划分为信号特征值(特征向量)与噪声特征值(特征向量)。

设λ_i是矩阵R_X的第i个特征值，v_i是与λ_i相对应的特征向量，则有

σ²v₁＝(AR_sA^H+σ²I)v_i (10)

将(10)式右边展开与左边比较，可得

A^Hv_i＝0 i＝K+1,K+2,...,M (11)

上式表明，噪声特征值所对应的特征向量与矩阵A的列向量正交，而A的各列与各信号源的方向相对应。因而利用来波信号的噪声特征向量就可以确定信号源的方向角。

步骤C、根据信号子空间和噪声子空间定义建立空间谱函数的表达式

用噪声特征向量为列，构造一个噪声矩阵E_n

E_n＝[v_K+1，v_K+2，...,v_M] (12)

获得空间谱函数P

该式中分母是信号向量和噪声矩阵的内积，当天线阵列的响应矢量a(θ_k)和噪声矩阵E_n的各列正交时，该分母为零，但由于随机噪声的存在，它实际上为一最小值，因此P有一尖峰，尖峰对应的角度即为信号的波达方向。

步骤三、基于求二阶导数的方法处理ADS-B信号的空间谱函数及空间谱函数谱峰搜索方法；

现有的MUSIC算法仅能找出步骤二中所得到的一个尖峰，若希望在两个或更多的ADS-B信号来波方向上能够实现对信号波达方向角的成功分辨，现有的MUSIC算法在特定的信噪比、采样率和阵列天线的孔径情况下，分辨率并不能达到要求的精度，不能满足实时获取信号波达方向角的需求。考虑到每个信号的来波方向在信号空间谱函数上体现的是一个波峰。在信号来波方向很近的情况下，空间谱函数上就只能出现一个波峰，而是实际上这个波峰是由与信源数量相同的多个波峰叠加组成的，因为采样间隔有限的关系，空间谱函数是离散函数，不能显示出多个波峰。由此发明人通过求空间谱函数的极值的方法来求信号的波达方向，空间谱函数的极值对应的方向角即为ADS-B信号的波达方向角。具体方法如下：

第k个信源的到达角θ_k对应的空间谱函数表示为p(θ_k)，因此空间谱函数P可以表示为

P＝[p(θ₁) p(θ₂) ... p(θ_k) ... p(θ_K)] (14)

由上式可知空间谱函数P是离散函数。由二阶导数定义可知，二阶导数的过零点即为函数的极值点。因此可以通过求离散函数二阶导数的方式来求空间谱函数P的极值。

按照导数定义，对于连续函数f(x)来说，一阶导数为

二阶导数为

其中，f_i'表示函数f(x)在x＝i处的导数；

对于离散函数u(x)来说，在x＝l处的一阶导数u'_l的表达式可以表示成如下形式

其中，表示在l点的邻域n内求和，且0<|x-l|≤n,x＝1,2,...,N,n＝1,2,...,N，N为自变量x的取值上限；

同理，由离散函数一阶导数定义可知离散函数u(x)在x＝l的二阶导数u"_l可以表示为

u'_l和u"_l分别是求导点x＝l的邻域内的2n+1个数据，n的取值直接影响到该离散函数的倒数的逼近程度。当2n+1个数值越接近时，导数越接近真实值，因此当n＝1时，u'_l和u"_l越接近真实的导数值，上式可以重写为

将上述求解结果应用于本发明的模型中，即基于MUSIC算法的改进方法对离散空间谱函数求二阶导数，得到二阶导数空间谱函数可以表示为

其中：

θ_k为空间谱函数P的自变量，表示到达角；

表示空间谱函数P在自变量为θ_k时对应的值；

N为θ_k的取值上限，取为360；

再对上述求得的二阶导数空间谱函数进行谱峰搜索求极大值，极大值对应的自变量的值即为ADS-B信号的来波方向角。

有益效果：

本发明基于线型阵列ADS-B信号传输模型，在传统MUSIC算法的基础上，提出了一种对空间谱函数求离散函数二阶导数，通过对二阶导数空间谱函数进行谱峰搜索来估算信号的来波方向角，在低信噪比、低快拍数情况下角度估计的分辨率相比于传统的MUSIC算法实现了提升，为ADS-B信号的人为干扰抑制提供了可行的方法。

具体的，

(1)本发明所建立的线型阵列ADS-B信号采集模型，为后续的空间谱估计步骤提供了方法。

(2)本发明所提出的基于改进MUSIC算法的高分辨率波达方向估计方法，能够提高传统的MUSIC算法的分辨精度，在低信噪比的情况下提高了传统MUSIC算法的分辨成功率。

(3)、本发明提出的谱峰搜索算法能够有效搜索到信号空间谱函数的谱峰。

附图说明

图1为本发明的流程图；

图2为线型ADS-B信号接收阵列示意图；

图3本传统MUSIC算法的流程图；

图4为传统MUSIC算法得到的接收空间谱函数中来波方向为30°和35°时功率谱示意；

图5为传统MUSIC算法得到的接收空间谱函数中来波方向为30°和31°时功率谱示意图；

图6为改进后的算法下的接收空间谱函数中来波方向为30°和35°时的功率谱示意图；

图7为改进后的算法下的接收空间谱函数中来波方向为30°和31°时的功率谱示意图；

图8为不同信噪比下算法成功分辨率；

图9为不同采样率下算法成功分辨率。

具体实施方式

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

下面结合图1至图7表明本实施方式，本实施方式包括以下步骤：

步骤一、按照表1中所列参数，模拟远场窄带信号参数和阵列天线。

表1

在此基础上，建立线型阵列的ADS-B信号采集模型，确定了线型阵列的阵列流型，设出信号矩阵和噪声矩阵，即可得到线型阵列的ADS-B信号采集模型，如图2所示。图2为线型ADS-B信号接收阵列示意图。

步骤二、估计ADS-B信号的空间谱函数：利用信号子空间和噪声子空间的正交性，先求出阵列信号的协方差矩阵，对协方差矩阵进行特征值分解，在对特征值进行排序，分辨出噪声子空间和信号子空间，最后得出空间谱函数。图3为ADS-B信号的空间谱函数估计流程图，图4和图5分别给出了传统MUSIC算法下得到的接收空间谱函数中不同来波方向的功率谱示意。

步骤三、基于求二阶导数的方法处理ADS-B信号的空间谱函数。如图6和图7为改进后的算法下得到的接收信号的不同波达方向时的空间谱函数示意图。

步骤四、利用蒙特卡洛仿真，对传统的MUSIC算法和本发明的基于求二阶导数的MUSIC算法分别求算法成功分辨的概率，设定成功分辨的判定条件为空间谱函数的极值的数量与波达方向的数量是否相等，如果相等则视为成功分辨；如果不相等，则视为不能成功分辨。

由前述理论可知算法的分辨率随着信号信噪比和采样快拍数的变化而变化，信噪比高，快拍数大，则算法的分辨率就高。但在实际中信号的信噪比和采样快拍数不能无限大，我们需要知道在小信噪比小快拍数的情况下算法的分辨性能。因此本例采用控制变量的方法，当仿真信号信噪比在一定范围内变化时采样快拍数不变，当采样快拍数在一定范围内变化的时候控制仿真信号信噪比不变。假设信号噪声为高斯白噪声，对两种方 法分别进行300次蒙特卡洛实验，将能分辨出到达角的估计次数统计求和，分别算出两种方法成功分辨的概率。

结果如图8和图9所示，图8为在不同条件下，两种算法成功分辨的概率，图8表示信噪比在一定氛围内变化的情况下两种算法成功分辨的概率，图9表示采样率在一定范围内变化的情况下两种算法成功分辨的概率从图上可以看出，在信噪比较小和采样率较小时，本发明的算法的成功分辨率明显高于传统算法，当信噪比达到12dB以上或者采样点数达到700次/秒以上时，两种方法分辨率相当且都接近100％。因此在信噪比小和采样率小的情况下，本发明具有明显优势。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (3)

1.一种基于改进的MUSIC算法的ADS-B信号波达方向估计方法，其特征在于：该方法包括顺序执行的以下步骤：

步骤一、对于模拟ADS-B信号进行简化，并对模拟天线阵列建立线型阵列的ADS-B信号采集模型；

步骤二、利用MUSIC算法估计ADS-B信号的空间谱函数P；

步骤三、根据获得的ADS-B信号的空间谱函数P，对空间谱函数P求解二阶导数获得二阶导数空间谱函数，并对获得的二阶导数空间谱函数进行谱峰搜索，谱峰对应的到达角即为ADS-B信号的波达方向。

2.根据权利要求1所述的一种基于改进的MUSIC算法的ADS-B波达方向估计方法，其特征在于：步骤一中，建立线型阵列的ADS-B信号采集模型的方法为：

设定接收信号为远场窄带信号，信源和天线阵列在同一平面内，入射到天线阵列的入射波为平面波；

以来波方向入射M根天线、阵元间距为d的均匀线阵，入射信号波长为λ，入射信号的信源总数为K个，第k个信源的到达角为θ_k，则天线阵列的响应矢量a(θ_k)为

$a\left({\mathrm{\&theta;}}_k\right)={\left[\begin{array}{cccc}1& \exp (-j2\mathrm{\&pi;}\frac{d}{\mathrm{\&lambda;}}{\mathrm{sin\&theta;}}_k)& \mathrm{...}& \exp (-j2\mathrm{\&pi;}(M-1\left)\frac{d}{\mathrm{\&lambda;}}{\mathrm{sin\&theta;}}_k\right)\end{array}\right]}^T$

则线型阵列的ADS-B信号采集模型表示为

X(t)＝AS(t)+N(t)

其中：

X(t)为阵元接收到的信号；

A为方向矩阵，且A＝[a(θ₁)，a(θ₂)，...,a(θ_K)]；

S(t)为信号向量；

N(t)为噪声向量。

3.根据权利要求1所述的一种基于改进的MUSIC算法的ADS-B波达方向估计方法，其特征在于：步骤三具体包括如下过程：

步骤A、对空间谱函数P按照下列计算式求解其二阶导数空间谱函数P"_θk

其中：

θ_k为空间谱函数P的自变量，表示到达角；

P_θk表示空间谱函数P在自变量为θ_k时对应的值；

N为θ_k的取值上限，取为360；

步骤B、对求得的二阶导数空间谱函数P"_θ进行谱峰搜索求极大值，极大值对应的角度即为ADS-B信号的波达方向。