基于经验模态分解的智能天线波达方向估计系统及方法

摘要

本发明公开了一种基于经验模态分解的智能天线波达方向估计系统及方法，系统包括天线阵模块，天线阵模块接收信号并传给经验模态分解模块，经验模态分解模块将信号分解成不同时间尺度的局部特征信号，并将局部特征信号中的高频信号传给高频方向估计模块，得到高频信号对应的阵列空间谱函数，将局部特征信号中的低频信号传给低频方向估计模块，得到低频信号对应的阵列空间谱函数，最后将高频信号对应的阵列空间谱函数和低频信号对应的阵列空间谱函数送至方向融合模块，由方向融合模块将两者对应累加，通过搜索累加后的阵列空间谱函数的峰值，确定信号的波达方向。本发明兼顾计算复杂度和分辨率，在复杂电磁环境下能够精确获取快信号目标。

说明书

技术领域

本发明涉及智能天线，特别是涉及一种基于经验模态分解的智能天线波达方向估计系统及方法。

背景技术

近年来，智能天线在提高系统通信质量、缓解频谱资源不足与无线通信日益发展的矛盾、以及降低系统造价和提高系统管理水平等方面，都具有独特的优点。目前智能天线技术不仅用于声纳、雷达、军事抗干扰通信，完成空间定位和滤波，还结合现代数字信号处理技术在基带形成天线波束，用于具有复杂电波传播环境的移动通信。经研究发现，智能天线可产生空间定向波束，使天线主瓣对准期望信号到达方向，零陷对准干扰信号到达方向，达到高效接收期望信号并抑制干扰信号的目的。

传统的智能天线波达方向估计方法有延迟-相加法(Bartlett)、多信号分类法(MUSIC)等。Bartlett法不需谱分解，计算简单，但分辨率较低。由于日益增长的使用需求和电磁环境的复杂性，现代电子测量装备需要比传统的探测设备具有更精准的信号目标侦测能力。因此分辨率较低的Bartlett法在实际应用中受到了限制。MUSIC法分辨率高、抗干扰效果好，但计算复杂度较高，实时效果差，不易捕获快信号目标。目前主要的研究思路有三种，第一种是寻找新的计算量小的波达方向估计方法，但这类方法的分辨率往往不及MUSIC算法；第二种是改进现有的波达方向估计方法，这类方法实际上是在计算量和精度之间进行了折中；第三种根据运算流程合理设计硬件，这类方法是以牺牲硬件成本为代价的，且效果有限。因此，需要一种计算复杂度和分辨率兼顾，满足复杂电磁环境下精确获取快信号目标的天线波达方向估计方法。

发明内容

发明目的：本发明的目的是提供一种兼顾计算复杂度和分辨率，在复杂电磁环境下能够精确获取快信号目标的基于经验模态分解的智能天线波达方向估计系统和方法。

技术方案：为达到此目的，本发明采用以下技术方案：

本发明所述的基于经验模态分解的智能天线波达方向估计系统，包括天线阵模块，天线阵模块接收信号并传给经验模态分解模块，经验模态分解模块将信号分解成不同时间尺度的局部特征信号，并将局部特征信号中的高频信号传给高频方向估计模块，得到高频信号对应的阵列空间谱函数，将局部特征信号中的低频信号传给低频方向估计模块，得到低频信号对应的阵列空间谱函数，最后将高频信号对应的阵列空间谱函数和低频信号对应的阵列空间谱函数送至方向融合模块，由方向融合模块将两者对应累加，通过搜索累加后的阵列空间谱函数的峰值，确定信号的波达方向。

本发明所述的基于经验模态分解的智能天线波达方向估计方法，包括以下步骤：

S1：天线阵模块接收信号并传给经验模态分解模块；

S2：经验模态分解模块将信号分解成不同时间尺度的局部特征信号，然后将局部特征信号中的高频信号传给高频方向估计模块，将局部特征信号中的低频信号传给低频方向估计模块；

S3：高频方向估计模块采用基于粒子群优化的多信号分类法得到高频信号对应的阵列空间谱函数；

S4：低频方向估计模块采用基于切比雪夫约束的延迟—相加法得到低频信号对应的阵列空间谱函数；

S5：将高频信号对应的阵列空间谱函数和低频信号对应的阵列空间谱函数送至方向融合模块，由方向融合模块将两者对应累加，通过搜索累加后的阵列空间谱函数的峰值，确定信号的波达方向。

进一步，所述步骤S3包括以下步骤：

S3.1：估计经验模态分解模块分解出的高频信号的协方差矩阵对协方差矩阵进行特征值分解，得到特征值和特征向量；

S3.2：确定协方差矩阵的最小特征值，得出最小特征值对应的特征向量，利用特征向量构造噪声特征向量矩阵E_N；

S3.3：随机产生种族规模为L的初始粒子群，根据式(1)更新每个粒子的空间谱函数：

$P_{MUSIC}\left({\theta }_i\right)=\frac{1}{a^H\left({\theta }_i\right)E_N{E}_N^{H}a\left({\theta }_i\right)}---\left(1\right)$

式(1)中，a(θ_i)为高频信号中第i个粒子对应的向量，P_MUSIC(θ_i)为第i个粒子的阵列空间谱函数，i＝1,2,...,L；

S3.4：计算每个粒子的阵列空间谱函数和目标空间谱函数之间的误差，根据式(2)设置误差最大值的最小化为粒子群的适应度函数fitness(x)：

fitness(x)＝min max(|P_MUSIC(θ_i)-P_d(θ_i)|²)>

式(2)中，P_d(θ_i)为第i个粒子的目标空间谱函数，当适应度函数fitness(x)达到最小值时的阵列空间谱函数记为P_G-MUSIC(θ)；

S3.5：对于n维空间，由L个粒子组成的种群，第i个粒子的速度V_i为：

V_i＝(v_i1,v_i2,...,v_in)(3)

式(3)中，v_id为第i个粒子速度在第d维空间的分量，i＝1,2,...,L，d＝1,2,...,n；

第i个粒子的位置X_i为：

X_i＝(x_i1,x_i2,...，x_in)>

式(4)中，x_id为第i个粒子位置在第d维空间的分量；

第i个粒子的个体极值P_i为：

P_i＝(p_i1,p_i2,...,p_in)(5)

式(5)中，p_id为第i个粒子个体极值在第d维空间的分量；

所有粒子的全局极值P_q为：

P_q＝(p_q1,p_q2,...,p_qn)(6)

式(6)中，p_qd为全局极值在第d维空间的分量；

根据式(7)和(8)更新粒子的位置及速度：

$v_{id}^{(t+1)}={\mathrm{wv}}_{id}^{\left(t\right)}+c_1{r}_1(p_{qd}^{\left(t\right)}-x_{id}^{\left(t\right)})+c_2{r}_2(p_{id}^{\left(t\right)}-x_{id}^{\left(t\right)})---\left(7\right)$

$x_{id}^{(t+1)}=x_{id}^{\left(t\right)}+v_{id}^{(t+1)}---\left(8\right)$

式(7)、(8)中，为t时刻第i个粒子速度在第d维空间的分量，为t+1时刻第i个粒子速度在第d维空间的分量，为t时刻第i个粒子位置在第d维空间的分量，为t+1时刻第i个粒子位置在第d维空间的分量，为t时刻第i个粒子个体极值在第d维空间的分量，为t时刻所有粒子的全局极值在第d维空间的分量，c₁、c₂为学习因子，r₁、r₂为(0,1)区间内的随机数，w为惯性权重；

S3.6：检查适应度函数fitness(x)是否达到最小值：如果达到最小值，则结束；否则，转至步骤S3.3继续执行。

进一步，所述步骤S3.5中惯性权重w的取值范围为0.3～0.8。

进一步，所述步骤S4包括以下步骤：

S4.1：根据式(9)将经验模态分解模块分解出的低频信号向量b(θ)点乘20dB的切比雪夫系数w_cheb，得到b_T：

b_T＝b(θ).×w_cheb(9)

S4.2：令W＝b_T，代入式(10)中得到式(11)：

P_out＝W^HR_XXW(10)

式(10)中，P_out为天线阵的总输出功率，R_XX为b(θ)的自相关函数；

然后，计算低频信号对应的阵列空间谱函数P_C-Bartlett(θ)：

P_C-Bartlett(θ)＝P_out＝(b(θ).×w_cheb)^HR_xx(b(θ).×w_cheb)>

进一步，所述步骤S5包括以下步骤：

S5.1：根据式(12)将高频信号对应的阵列空间谱函数和低频信号对应的阵列空间谱函数对应累加，得到方向融合后的阵列空间谱函数P_EMD-GMCB(θ)；

P_EMD-GMCB(θ)＝P_G-MUSIC(θ)+P_C-Bartlett(θ)(12)；

S5.2：通过搜索方向融合后的阵列空间谱函数P_EMD-GMCB(θ)的峰值，确定信号波达方向。

有益效果：

1)本发明克服了延迟—相加法(Bartlett)分辨率低、副瓣较多的缺点，副瓣幅度平均降低9dB，约降低总幅度的62.9％，较Bartlett法运行时间减少了37.5％，从而提高了系统对干扰信号的抑制力以及获取期望信号的准确性；

2)本发明克服了多信号分类法(MUSIC)计算复杂度高的缺点，降低MUSIC法95.2％的计算复杂度，较MUSIC法运行时间短，平均运行时间约减少了35.5％；

3)本发明实现了分辨率和计算复杂度的高度兼顾，满足复杂性电磁环境里精确获取快信号目标的使用需求。

附图说明

图1为本发明的系统框图；

图2为本发明的具体结构图；

图3为本发明的具体实施方式中天线阵模块接收到的信号的波形图；

图4为本发明的具体实施方式中经验模态分解模块分解得到的局部特征信号中的第一高频信号imf1的波形图；

图5为本发明的具体实施方式中经验模态分解模块分解得到的局部特征信号中的第二高频信号imf2的波形图；

图6为本发明的具体实施方式中经验模态分解模块分解得到的局部特征信号中的第三高频信号imf3的波形图；

图7为本发明的具体实施方式中经验模态分解模块分解得到的局部特征信号中的第四高频信号imf4的波形图；

图8为本发明的具体实施方式中经验模态分解模块分解得到的局部特征信号中的第一低频信号imf5的波形图；

图9为本发明的具体实施方式中经验模态分解模块分解得到的局部特征信号中的第二低频信号imf6的波形图；

图10为本发明的具体实施方式中经验模态分解模块分解得到的局部特征信号中的第三低频信号imf7的波形图；

图11为本发明的具体实施方式中经验模态分解模块分解得到的局部特征信号中的第四低频信号imf8的波形图；

图12为本发明的具体实施方式中采用的基于粒子群优化的多信号分类法和传统的多信号分类法得到的高频信号的阵列空间谱函数的对比图；

图13为本发明的具体实施方式中采用的基于切比雪夫约束的延迟—相加法和传统的延迟—相加法得到的低频信号的阵列空间谱函数的对比图；

图14为本发明的具体实施方式得到的累加后的阵列空间谱函数与分别采用传统的多信号分类法以及延迟—相加法得到的阵列空间谱函数的对比图。

具体实施方式

下面结合具体实施方式对本发明的技术方案作进一步的介绍。

本发明公开了一种基于经验模态分解的智能天线波达方向估计系统，如图1和图2所示，包括天线阵模块1，天线阵模块1接收信号并传给经验模态分解模块2，经验模态分解模块2将信号分解成不同时间尺度的局部特征信号，并将局部特征信号中的高频信号传给高频方向估计模块3，得到高频信号对应的阵列空间谱函数，将局部特征信号中的低频信号传给低频方向估计模块4，得到低频信号对应的阵列空间谱函数，最后将高频信号对应的阵列空间谱函数和低频信号对应的阵列空间谱函数送至方向融合模块5，由方向融合模块5将两者对应累加，通过搜索累加后的阵列空间谱函数的峰值，确定信号的波达方向。

本发明还公开了一种基于经验模态分解的智能天线波达方向估计方法，步骤如下：

S1：天线阵模块1接收信号并传给经验模态分解模块2；假设期望信号s_k(t)从Q_k方向入射，干扰信号s_j(t)从Q_j方向入射，其中干扰信号共计k-1个，则天线阵元接收的信号x(t)可表示为：

$x\left(t\right)=S_k\left(t\right)+\underset{j=1}{\overset{M-1}{\Sigma }}{S}_j\left(t\right)+n\left(t\right)---\left(1\right)$

式(1)中，n(t)为噪声信号，S_k(t)和S_j(t)分别为：

S_k(t)＝s_k(t)·V(Q_k)S_j(t)＝s_j(t)·V(Q_j)>

S2：经验模态分解模块2将信号分解成不同时间尺度的局部特征信号，也即在matlab中，采用EMD函数，用语句imf＝emd(X)，将阵列信号X分解到多时间尺度频域；用语句plot(c,X)；title('rignal')，可得到原始信号图形，用语句subplot(4,1,2)；plot(c,imf{1})；title('imf1')，可得到局部特征信号imf1，同理可得局部特征信号imf2、imf3、imf4、imf5、imf6、imf7、imf8，其中imf1、imf2、imf3、imf4为高频信号，如图4—图7所示，imf5、imf6、imf7、imf8为低频信号，如图8—图11所示；从而能够实现信号的正交处理，可减小自适应滤波器的输入向量自相关阵特征值的分散程度，大大增加算法的收敛步长，提高收敛速度；此外，经验模态分解模块2还将高频信号imf1、imf2、imf3和imf4传给高频方向估计模块3，将低频信号imf5、imf6、imf7、imf8传给低频方向估计模块4；

S3：高频方向估计模块3采用基于粒子群优化的多信号分类法得到高频信号imf1、imf2、imf3、imf4对应的阵列空间谱函数，其中，基于粒子群优化的多信号分类法通过L个粒子组成的粒子群优化多信号分类法得到，包括以下步骤：

S3.1：估计经验模态分解模块2分解出的高频信号的协方差矩阵对协方差矩阵进行特征值分解，得到特征值和特征向量；

S3.2：确定协方差矩阵的最小特征值，得出最小特征值对应的特征向量，利用特征向量构造噪声特征向量矩阵E_N；

S3.3：随机产生种族规模为L的初始粒子群，根据式(3)更新每个粒子的空间谱函数：

$P_{MUSIC}\left({\theta }_i\right)=\frac{1}{a^H\left({\theta }_i\right)E_N{E}_N^{H}a\left({\theta }_i\right)}---\left(3\right)$

式(3)中，a(θ_i)为高频信号中第i个粒子对应的向量，P_MUSIC(θ_i)为第i个粒子的阵列空间谱函数，i＝1,2,...,L；

S3.4：计算每个粒子的阵列空间谱函数和目标空间谱函数之间的误差，根据式(4)设置误差最大值的最小化为粒子群的适应度函数fitness(x)：

fitness(x)＝min max(|P_MUSIC(θ_i)-P_d(θ_i)|²)>

式(4)中，P_d(θ_i)为第i个粒子的目标空间谱函数，当适应度函数fitness(x)达到最小值时的阵列空间谱函数记为P_G-MUSIC(θ)；

S3.5：对于n维空间，由L个粒子组成的种群，第i个粒子的速度V_i为：

V_i＝(v_i1,v_i2,...,v_in)>

式(5)中，v_id为第i个粒子速度在第d维空间的分量，i＝1,2,...,L，d＝1,2,...,n；

第i个粒子的位置X_i为：

X_i＝(x_i1,x_i2,...，x_in)(6)

式(6)中，x_id为第i个粒子位置在第d维空间的分量；

第i个粒子的个体极值P_i为：

P_i＝(p_i1,p_i2,...,p_in)(7)

式(7)中，p_id为第i个粒子个体极值在第d维空间的分量；

所有粒子的全局极值P_q为：

P_q＝(p_q1,p_q2,...,p_qn)>qd为全局极值在第d维空间的分量；

根据式(9)和(10)更新粒子的位置及速度：

$v_{id}^{(t+1)}={\mathrm{wv}}_{id}^{\left(t\right)}+c_1{r}_1(p_{qd}^{\left(t\right)}-x_{id}^{\left(t\right)})+c_2{r}_2(p_{id}^{\left(t\right)}-x_{id}^{\left(t\right)})---\left(9\right)$

$x_{id}^{(t+1)}=x_{id}^{\left(t\right)}+v_{id}^{(t+1)}---\left(10\right)$

式(9)、(10)中，为t时刻第i个粒子速度在第d维空间的分量，为t+1时刻第i个粒子速度在第d维空间的分量，为t时刻第i个粒子位置在第d维空间的分量，为t+1时刻第i个粒子位置在第d维空间的分量，为t时刻第i个粒子个体极值在第d维空间的分量，为t时刻所有粒子的全局极值在第d维空间的分量，c₁、c₂为学习因子，r₁、r₂为(0,1)区间内的随机数，w为惯性权重，w的取值范围为0.3～0.8；

S3.6：检查适应度函数fitness(x)是否达到最小值：如果达到最小值，则结束；否则，转至步骤S3.3继续执行；

S4：低频方向估计模块4采用基于切比雪夫约束的延迟—相加法得到低频信号imf5、imf6、imf7、imf8对应的阵列空间谱函数，其中，基于切比雪夫约束的延迟—相加法通过切比雪夫系数约束延迟—相加法，包括以下步骤：

S4.1：根据式(11)将经验模态分解模块2分解出的低频信号向量b(θ)点乘20dB的切比雪夫系数w_cheb，得到b_T：

b_T＝b(θ).×w_cheb(11)

S4.2：令W＝b_T，代入式(12)中得到式(13)：

P_out＝W^HR_XXW(12)式(12)中，P_out为天线阵的总输出功率，R_XX为b(θ)的自相关函数；

然后，计算低频信号对应的阵列空间谱函数P_C-Bartlett(θ)：

P_C-Bartlett(θ)＝P_out＝(b(θ).×w_cheb)^HR_xx(b(θ).×w_cheb)(13)。

S5：将高频信号imf1、imf2、imf3、imf4对应的阵列空间谱函数和低频信号imf5、imf6、imf7、imf8对应的阵列空间谱函数送至方向融合模块5，由方向融合模块5将两者对应累加，通过搜索累加后的阵列空间谱函数的峰值，确定信号的波达方向，具体步骤如下：

S5.1：根据式(14)将高频信号imf1、imf2、imf3、imf4对应的阵列空间谱函数和低频信号imf5、imf6、imf7、imf8对应的阵列空间谱函数对应累加，得到方向融合后的阵列空间谱函数P_EMD-GMCB(θ)；

P_EMD-GMCB(θ)＝P_G-MUSIC(θ)+P_C-Bartlett(θ)>

S5.2：通过搜索方向融合后的阵列空间谱函数P_EMD-GMCB(θ)的峰值，确定信号波达方向。

本具体实施方式中，天线阵模块1采用天线阵元数目为20的直线阵，天线阵模块1接收到入射方向分别为-10°、20°的两个均值为零的信号，如图3所示。

本具体实施方式对多信号分类法(MUSIC法)和基于粒子群优化的多信号分类法(G-MUSIC法)进行了比较，分别采用这两种方法得到了高频信号的阵列空间谱函数，如图12所示。具体比较过程如下：

A)用Matlab仿真MUSIC法实验。设入射方向分别为-10°，20°，采样数为K＝500个快拍，信噪比为SNR＝10dB，角度扫描范围是[0°,180°]，扫描间隔为1°。采用多信号分类法(MUSIC)进行阵列空间谱函数搜索，阵列空间谱函数如图12所示。

B)用Matlab仿真G-MUSIC法实验。设粒子数L为20，粒子最大速率为vm max1为40，vm max1为190，最大惯性因子为0.8，最小惯性因子为0.1，学习因子c1、c2都为2，最大进化代数为300。设入射方向分别为-10°，20°，采样数设为K＝500个快拍，信噪比设为SNR＝10dB，角度扫描范围是[0°,180°]，扫描间隔为1°。采用粒子群优化多信号分类法(G-MUSIC)进行阵列空间谱函数搜索，阵列空间谱函数如图12所示。

C)对比分析MUSIC法与G-MUSIC法。从图12中可以看出，通过搜索阵列空间谱函数的峰值，MUSIC法与G-MUSIC法均可确定信号波达方向分别为-10°、20°。两种方法的峰值尖锐且都没有旁瓣，算法分辨率都较高，G-MUSIC法较MUSIC法主瓣宽一点。多信号分类法(MUSIC)，其整个空间搜索需要适应度计算的次数为360×90＝32400次，而基于粒子群的多信号分类法(G-MUSIC)，其整个空间搜索需要适应度计算的次数为80×20＝1600次。因此G-MUSIC法的计算复杂度远远小于MUSIC法的计算复杂度，G-MUSIC法仅仅约有MUSIC计算量的5％，计算复杂度得到了较大程度的改善。

本具体实施方式对延迟—相加法(Bartlett法)和基于切比雪夫约束的延迟—相加法(C-Bartlett法)进行了比较，分别采用这两种方法得到了低频信号的阵列空间谱函数，如图13所示。具体比较过程如下：

A)用Matlab仿真Bartlett法实验。设入射方向分别为-10°，20°，采样数为K＝500个快拍，信噪比为SNR＝10dB。角度扫描范围是[0°,180°]，扫描间隔为1°。采用延迟—相加法(Bartlett)进行阵列空间谱函数搜索，阵列空间谱函数如图13所示。

B)用Matlab仿真C-Bartlett法实验。两信源入射方向分别为-10°，20°，采样数设为K＝500个快拍，信噪比设为SNR＝10dB。角度扫描范围是[0°,180°]，扫描间隔为1°。设切比雪夫系数w_cheb为-20dB，采用切比雪夫约束延迟—相加法(C-Bartlett)进行阵列空间谱函数搜索，阵列空间谱函数如图13所示。

C)对比分析Bartlett法与C-Bartlett法。由图13可以看出，Bartlett法中随着角度的改变输出功率有两个峰值，分别对应两个信号的入射方向-10°，20°，但是峰值不够尖锐且旁瓣较多，可见Bartlett算法分辨率不是很高。但经过切比雪夫约束后的C-Bartlett法实现了入射信号的方向估计，且副瓣个数减少，副瓣幅度平均降低9dB，约降低总幅度的62.9％，主瓣也比普通的Bartlett法稍尖锐，分辨率得到改善。

此外，本具体实施方式还将本发明方法(EMD-GMCB法)与传统的MUSIC法以及Bartlett法进行了比较，结果如图14所示。具体比较过程如下：

A)计算量对比

记MUSIC法、Bartlett法计算复杂度分别为O₁(M)、O₂(M)，其中O₁(M)＞＞O₂(M)。记G-MUSIC法、C-Bartlett法计算复杂度分别为O′₁(M)、O'₂(M)，资料表明O′₁(M)＞O′₂(M)。前文表明G-MUSIC法仅仅约有MUSIC计算量的5％，记为O′₁(M)≈5％O₁(M)。又很据imf1、imf2、imf3、imf4、imf5、imf6、imf7、imf8能量大致分布，可知较高频部分信号约占总能量的95％，较低频信号约占总能量的5％。则可得EMD-GMCB法在一次迭代中的计算复杂度约为95％×O′₁(M)+5％×O'₂(M)≈95％×O′₁(M)≈95％×5％O₁(M)≈4.8％O₁(M)，记为O₃(M)。因此可得O₃(M)＜＜O₁(M)，EMD-GMCB法的计算复杂度O₃(M)远远低于MUSIC法的复杂度O₁(M)，降低MUSIC法95.2％的计算复杂度，可见本发明提出的EMD-GMCB法是一种计算复杂度较低的方法。

B)分辨率对比

相同条件下，Bartlett、MUSIC、EMD-GMCB法空间谱图如图14所示。可见，EMD-GMCB法克服Bartlett法副瓣较高的缺点，副瓣幅度平均降低9dB，约降低总幅度的62.9％，从而提高了系统对干扰信号的抑制力以及获取期望信号的准确性。根据计算量分析，EMD-GMCB法计算复杂度为O₃(M)≈4.8％O₁(M)，即O₃(M)＜＜O₁(M)，降低MUSIC法95.2％的计算复杂度。可见本发明提出的EMD-GMCB法实现了分辨率和计算复杂度的高度兼顾。

C)运行时间对比

本实验在联想电脑G460AP机型运行Bartlett法、MUSIC法、EMD-GMCB法相应的Matlab程序，三种方法运行时间如表1所示。由表1可看出，在相同阵元数目M下，EMD-GMCB法运行时间最短，较Bartlett法运行时间提高了37.5％，较MUSIC法运行时间提高了35.5％。运行时间既能反应计算复杂度也能反应收敛性，可见本发明提出的EMD-GMCB法是一种收敛性较高的方法。

表1Bartlett算法、MUSIC法、EMD-GMCB法运行时间