CN108520195B - 一种基于引力搜索算法的music谱峰搜索方法 - Google Patents

Abstract

本发明提出一种基于引力搜索算法的MUSIC谱峰搜索方法，该方法将GSA算法和MUSIC算法相结合，具体方法是将MUSIC算法中的空间谱函数设为GSA算法的适应度函数，利用GSA算法优秀的全局寻优能力，迅速地找到空间谱函数的极值点，完成了对谱函数的空间搜索。仿真表明，该方法取得了较好的估计效果，在实验条件相同的情况下，与MUSIC算法中常用的均匀遍历法相比，该方法使得计算量大大降低。

Description

一种基于引力搜索算法的MUSIC谱峰搜索方法

技术领域

本发明涉及电子通信技术领域，具体涉及一种基于引力搜索算法的MUSIC谱峰搜索方法。

背景技术

阵列信号处理是电子通信领域的重点方向之一，在雷达、声纳和语音信号处理相关领域得到了广泛应用，最近以来引起了很高的重视，其发展重心主要集中波束形成和DOA估计的问题上。阵列信号处理是指在不同位置的多个传感器按特定的规则组成传感器阵列，利用传感器阵列估计信号的来波方向并且通过波束形成的方法对接收的信号进行处理，随着高分辨率谱估计技术的发展，DOA估计的问题也愈加受到重视。

经典的DOA估计方法中，以MUSIC(Multiple Signal Classification，多重信号分类)算法的应用最为广泛，该算法具有高精度和稳定性好的优点。MUSIC算法是一种典型的子空间结构类算法，其首先将接收数据构造的协方差矩阵分解为噪声子空间和信号子空间，然后利用噪声子空间与信号子空间的方向矢量之间的相互正交性来构造空间谱函数，谱函数对应的极值点即为参数的估计值。由于MUSIC算法具有较高的估计精度和稳定性，且突破了瑞利限制，目前已成为最具研究价值的DOA算法之一。传统的MUSIC算法中通常采用均匀遍历法计算谱函数，即对整个搜索空间进行遍历搜索，估计精度取决于间隔点间的步长大小，虽然该方法能将所有方向都能搜索到，但当精度要求较高时计算量也会相应变大，特别是对于多维、多目标的参数估计情况，MUSIC算法的谱峰搜索更加复杂，所需的时间也更多。由于其谱峰搜索部分的计算复杂度较高，使得该算法在实时性要求较高的场合的应用中受到限制。

MUSIC算法中谱峰搜索部分的计算量往往占该算法总计算量的一半以上，为了提高MUSIC算法的运算速度，使其满足实时性的要求，就要选取合适的方法进行MUSIC谱峰搜索。本发明主要针对MUSIC算法的谱峰搜索方法进行创新，在满足精度要求情况下提高搜索的速度。

下面对信号模型、MUSIC算法原理、引力搜索算法原理进行简单介绍。

1信号模型

假设阵列的信号接收模型为均匀线阵，且阵列天线由M个传感器构成，相邻传感器间隔为d，阵列外有D个波长为λ的窄带平面波，从不同入射角度入射，相邻传感器间隔d＝λ/2。则阵列接收数据的矢量形式为：

X(t)＝A(θ)S(t)+n(t) (1)

其中，X(t)＝[x₁(t),x₂(t),...,x_M(t)]^T为M维数据接收向量，S(t)＝[s₁(t),s₂(t),...,s_D(t)]^T为D维的信号源向量，n(t)是空间和时间均独立的高斯白噪声矢量，均值为0，方差为σ²。阵列流型矩阵A(θ)为M×D阶矩阵，其表达式为：

A(θ)＝[a(θ₁),a(θ₂),...a(θ)_D] (2)

其中，第k列a(θ_k)为第k个信号源s_k的导向向量，即：

其中，j指虚数单位，

θ_k为第k个信号源的来波方向角。阵列数据的协方差矩阵为：

R＝E[X(t)X^H(t)]＝AR_SA+R_N (4)

式中，E[·]和[·]^H分别指期望值和共轭矩阵，A和X(t)分别指阵列流型矩阵和数据接收向量，R_S和R_N分别为信源的协方差矩阵和噪源的协方差矩阵。

2MUSIC算法原理

MUSIC算法是基于噪声子空间类的算法，它利用信号子空间和噪声子空间的正交性，构造空间谱函数，通过谱峰搜索，估计信号的来波方向。MUSIC算法是基于以下假设的数学模型上进行的：

假设1：阵列类型为各向同性均匀线阵，且阵元间隔一般小于最高频率信号波长的二分之一；

假设2：空间中的噪声为高斯白噪声，不同阵元间的噪声均为平稳随机过程，且相互独立，各阵元间的噪声均值为零、方差相等；

假设3：空间信号为窄带平面波，它与阵元噪声相互独立；

假设4：信号源数小于阵元数。

在此假设基础上，MUSIC算法对波达方向DOA的估计在理论上可以有任意高的分辨率。由于MUSIC算法的超分辨特点，该算法及其改进算法在近几十年来受到了广泛关注，并已成为超分辨子空间算法中的经典算法之一。

MUSIC算法的实现步骤为：

(1)根据传感器阵列接收的N次快拍数据构造数据协方差矩阵的估计值如下：

式中X(i)为第i次接收到的数据。

(2)对数据协方差矩阵进行特征值分解：

R＝VΛV^H (6)

其中，V是与R同阶的方阵，Λ是由特征值组成的对角矩阵。

(3)对分解得到的特征值按大小进行排序，把与信号个数D相等的特征值所对应的特征向量作为信号子空间，把剩下的(M-D)个特征值对应的特征向量作为噪声子空间：R＝V_SΛ_SV_S ^H+V_NΛ_NV_N ^H，V_S和V_N为别是信号子空间和噪声子空间，并确定R的最小特征值个数n_E。

(4)构造空间谱函数：

其中，a(θ)为信号子空间的方向向量。

对谱函数进行空间搜索，其D个极大值对应的θ就是待估信源方向的估计值。

3引力搜索算法介绍

引力搜索算法(Gravitational Search Algorithm，GSA)是2009年由E.Rashedi根据万有引力现象和牛顿力学定理总结得出的一个仿生智能算法，该算法思想简单，搜索速度快，是全局优化问题的有效解决方案之一。

引力搜索算法(GSA)是一个较为常见的现代元启发式算法，根据算法的描述，有质量的个体之间产生相互作用力，质量较小的个体被吸引并产生加速度移动。每个个体有四个重要的参数：位置、惯性质量、主动惯性质量和被动惯性质量。惯性质量、主动惯性质量和被动惯性质量在算法中被认为是相等的。

与其他的群智能优化算法类似，该算法首先初始化种群的粒子，各个粒子之间存在相互作用力，随着算法的进行，这种作用力会促使所有质量小的粒子朝着相对而言质量更大的粒子移动。各粒子的速度和位置从初代开始改变，粒子所在的位置就是问题的可能解，每个粒子的性能与其惯性质量成一定的比例，当算法终止时，选择性能最好的粒子作为问题的解。

发明内容

本发明的目的是：提出一种基于引力搜索算法的MUSIC谱峰搜索方法，以期能降低运算复杂度，加快算法的执行速度，从而更好地应用于实际场景。

为了达到上述目的，本发明所采用的技术方案是：一种基于引力搜索算法的MUSIC谱峰搜索方法，其特征在于，将引力搜索算法用于MUSIC算法的谱峰搜索问题，则种群中粒子的位置代表着估计方位角在每次迭代过程中的解，所述方法包括如下步骤：

步骤一，初始化所有粒子在给定维度的速度和位置信息，并设置迭代次数和引力搜索算法中的初始参数，将MUSIC空间谱函数P_MUSIC设为引力搜索算法的适应度函数，将粒子的位置作为自变量θ，代入MUSIC空间谱函数P_MUSIC算出所有粒子的适应度值，并根据适应度值衡量粒子的性能；

步骤二，运用引力搜索算法在每一次迭代过程中更新重力常数，并利用适应度值算出每个粒子的惯性质量，根据重力常数和惯性质量计算出每个粒子的引力之和，再由引力之和和惯性质量得出每个粒子的加速度；

步骤三，利用加速度计算出每次迭代过程中所有粒子的速度和位置，适应度值最大的粒子的位置代表着DOA估计问题的解；

步骤四，检查算法的终止条件，如果满足，终止算法并选择适应度值最大的一组解作为DOA估计问题的解。

进一步地，所述基于引力搜索算法的MUSIC谱峰搜索方法的具体实现过程包括：

步骤1：随机地初始化种群中的粒子，每个粒子都有自己的初始位置，各粒子的位置X_i的表达式为：

其中，

和D分别表示粒子i在第k维搜索空间的位置和搜索空间的维数，D也等于待估计信源的个数，N为初始种群大小，算法在各搜索空间的搜索范围在[-90°,90°]之间；

步骤2：算出各粒子在不同时刻的惯性质量，计算公式如下：

M_ai＝M_pi＝M_i,i＝1,2,..,N (9)

M_ai和M_pi分别代表第i个粒子主动惯性质量和被动惯性质量，M_i(t)代表第i个粒子的惯性质量，fit_i(t)代表t时刻第i个粒子的适应度值，是根据MUSIC空间谱函数

得出的值，其中，a(θ)为信号子空间的方向向量，V_N为噪声子空间，[·]^H指共轭矩阵；

依据上述MUSIC空间谱函数P_MUSIC公式计算每次迭代中所有粒子的适应度，由于该问题是优化最大值问题，时刻t的最优解best(t)和最差解worst(t)分别为：

步骤3：从引力定律中可以推出，第k维上一个粒子i在t时刻受到另一个粒子j的作用力的计算公式如下：

R_ij(t)表示在t时刻两个粒子i和j的欧几里得度量，ε表示一个非常小的常数；

在t时刻的万有引力常数G(t)为：

G(t)＝G₀×e^-α·t/T (18)

G₀表示在t₀时刻G的取值，一般取值为100，α为下降系数，一般取值为20，t代表当前循环次数，T代表最大循环次数；

则在t时刻，X_i在第k维上受到的引力总和是由种群中其余粒子对其产生的全部引力相加得到，计算公式如下所示：

式中，rand_j表示在[0,1]之间的随机数；

步骤4：从牛顿的力学原理可以得知，粒子产生的加速度等于该粒子受到的引力总和与其惯性质量之比，则粒子i的加速度的计算公式如下：

然后在每一次更新过程中，根据计算得到的加速度来更新每个粒子的速度，再根据粒子的速度更新每个粒子的位置，计算公式如下所示：

与现有技术相比，本发明的有益效果是：本发明提出的基于引力搜索算法的MUSIC谱峰搜索方法，该方法将GSA算法和MUSIC算法相结合，具体方法是将MUSIC算法中的空间谱函数设为GSA算法的适应度函数，利用GSA算法优秀的全局寻优能力，迅速地找到空间谱函数的极值点，完成了对谱函数的空间搜索。仿真表明，该方法取得了较好的估计效果，在实验条件相同的情况下，与MUSIC算法中常用的均匀遍历法相比，该方法使得计算量大大降低。

附图说明

图1为本发明基于引力搜索算法的MUSIC谱峰搜索方法的流程图。

图2为本发明实施例的迭代次数与DOA估计角度的变化曲线。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

基于引力搜索算法的MUSIC谱峰搜索步骤：

将GSA算法用于MUSIC算法的谱峰搜索问题，则种群中粒子的位置代表着估计方位角在每次迭代过程中的解，该方法的具体描述如下：

步骤1：随机地初始化种群中的粒子，每个粒子都有自己的初始位置，各粒子的位置X_i的表达式为：

其中，

和D分别表示粒子i在第k维搜索空间的位置和搜索空间的维数，D也等于待估计信源的个数，N为初始种群大小。算法在各搜索空间的搜索范围在[-90°,90°]之间。

步骤2：算出各粒子在不同时刻的惯性质量，计算公式如下：

M_ai＝M_pi＝M_i,i＝1,2,..,N (9)

M_ai和M_pi分别代表第i个粒子主动惯性质量和被动惯性质量，M_i(t)代表第i个粒子的惯性质量，fit_i(t)代表t时刻第i个粒子的适应度值，是根据MUSIC空间谱函数

得出的值。

依据公式(7)计算每次迭代中所有粒子的适应度，由于该问题是优化最大值问题，时刻t的最优解best(t)和最差解worst(t)分别为：

步骤3：从引力定律中可以推出，第k维上一个粒子i在t时刻受到另一个粒子j的作用力

的计算公式如下：

R_ij(t)表示在t时刻两个粒子i和j的欧几里得度量，ε表示一个非常小的常数，

表示粒子i在第k维度t时刻的位置。在t时刻的万有引力常数G(t)为：

G(t)＝G₀×e^-α·t/T (15)

G₀表示在t₀时刻G的取值，一般取值为100，α为下降系数，一般取值为20，t代表当前循环次数，T代表最大循环次数。

则在t时刻，X_i在第k维上受到的引力总和F_i ^k(t)是由种群中其余粒子对其产生的全部引力相加得到，计算公式如下所示：

式中，

为在第k维度中粒子j对粒子i的作用力，rand_j表示在[0,1]之间的随机数。

步骤4：从牛顿的力学原理可以得知，粒子产生的加速度等于该粒子受到的引力总和与其惯性质量之比，则第k维度的粒子i在t时刻的加速度

的计算公式如下：

然后在每一次更新过程中，根据计算得到的加速度来更新每个粒子的速度，再根据粒子的速度更新每个粒子的位置，则第k维度中粒子i的速度与位置更新公式如下：

式中

和

分别表示第i个粒子在第k维度的位置和速度。

下面给出该方法的总体思路：

首先初始化所有粒子在各维度的速度和位置信息，并设置迭代次数与算法中的其它参数，MUSIC空间谱函数

设为GSA算法的适应度函数。将粒子的位置作为自变量θ，并代入P_MUSIC算出所有粒子的适应度值。适应度值是衡量粒子性能的重要指标。

然后在每一次迭代过程中更新重力常数，并利用适应度值算出每个粒子的惯性质量，根据重力常数和惯性质量求出每个粒子的引力之和，再由引力之和和惯性质量得到每个粒子的加速度。

最后利用加速度计算出每次迭代过程中所有粒子的速度和位置，由于本问题是优化最大值问题，所以适应度值最大的粒子的位置代表着DOA估计问题的解。

检查算法的终止条件，如果满足，终止算法并选择公式(19)中适应度值最大的一组解作为DOA估计问题的解。整个方法的流程图如图1所示。

仿真结果

仿真一：假设阵列传感器使用10元等距线阵，相邻传感器的间隔为来波信号波长的一半。并考虑空间中存在一个窄带平面波，其入射方向角是30^°。信号的快拍数为100，信噪比(Signal-Noise Ratio,SNR)固定为10dB。在GSA算法中，各参数设置如下：初始种群大小N＝10，G₀＝100，α＝20，最大循环次数T＝50。

图2给出了迭代次数与DOA估计角度的变化曲线。

由图2可以看到，在参数为一维的情况下，GSA算法迭代到第6次时即停止并成功估计出了单个来波信号的DOA，初始种群大小设置为10，则适应度函数公式(7)计算了10×6＝60次，相比较传统谱峰搜索方法的1×180＝180次，该方法使得计算量大大降低。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (2)

1.一种基于引力搜索算法的MUSIC谱峰搜索方法，应用于电子通信领域的阵列信号处理的DOA估计的问题上，解决MUSIC算法谱峰搜索部分的计算复杂度高的问题，其特征在于：

首先，定义信号模型：

假设阵列的信号接收模型为均匀线阵，且阵列天线由M个传感器构成，相邻传感器间隔为d，阵列外有D个波长为λ的窄带平面波，从不同入射角度入射，D表示信号个数，相邻传感器间隔d＝λ/2，则阵列接收数据的矢量形式为：

X(t)＝A(θ)S(t)+n(t) (1)

其中，t为时刻，X(t)＝[x₁(t),x₂(t),...,x_M(t)]^T为M维数据接收向量，S(t)＝[s₁(t),s₂(t),...,s_D(t)]^T为D维的信号源向量，n(t)是空间和时间均独立的高斯白噪声矢量，均值为0，方差为σ²，阵列流型矩阵A(θ)为M×D阶矩阵，其表达式为：

A(θ)＝[a(θ₁),a(θ₂),...a(θ)_D] (2)

其中，第k列a(θ_k)为第k个信号源s_k的导向向量，即：

其中，j指虚数单位，

θ_k为第k个信号源的来波方向角，阵列数据的协方差矩阵为：

R＝E[X(t)X^H(t)]＝AR_SA+R_N (4)

式中，E[·]和[·]^H分别指期望值和共轭矩阵，A和X(t)分别指阵列流型矩阵和数据接收向量，R_S和R_N分别为信源的协方差矩阵和噪源的协方差矩阵；

其次，关于MUSIC算法：

MUSIC算法是基于噪声子空间类的算法，它利用信号子空间和噪声子空间的正交性，构造空间谱函数，通过谱峰搜索，估计信号的来波方向，MUSIC算法是基于以下假设的数学模型上进行的：

假设1：阵列类型为各向同性均匀线阵，且阵元间隔一般小于最高频率信号波长的二分之一；

假设2：空间中的噪声为高斯白噪声，不同阵元间的噪声均为平稳随机过程，且相互独立，各阵元间的噪声均值为零、方差相等；

假设3：空间信号为窄带平面波，它与阵元噪声相互独立；

假设4：信号源数小于阵元数；

在此假设基础上，MUSIC算法的实现步骤为：

(1)根据传感器阵列接收的N次快拍数据构造数据协方差矩阵的估计值如下：

式中X(i)为第i次接收到的数据；

(2)对数据协方差矩阵进行特征值分解：

R＝VΛV^H (6)

其中，V是与R同阶的方阵，Λ是由特征值组成的对角矩阵；

(3)对分解得到的特征值按大小进行排序，把与信号个数D相等的特征值所对应的特征向量作为信号子空间，把剩下的(M-D)个特征值对应的特征向量作为噪声子空间：R＝V_SΛ_SV_S ^H+V_NΛ_NV_N ^H，V_S和V_N为别是信号子空间和噪声子空间，并确定R的最小特征值个数n_E；

(4)构造空间谱函数：

其中，a(θ)为信号子空间的方向向量；

对谱函数进行空间搜索，其D个极大值对应的θ就是待估信源方向的估计值；

再次，关于引力搜索算法：

引力搜索算法是一个现代元启发式算法，根据算法的描述，有质量的个体之间产生相互作用力，质量较小的个体被吸引并产生加速度移动；每个个体有四个重要的参数：位置、惯性质量、主动惯性质量和被动惯性质量，惯性质量、主动惯性质量和被动惯性质量在算法中被认为是相等的；

该算法首先初始化种群的粒子，各个粒子之间存在相互作用力，随着算法的进行，这种作用力会促使所有质量小的粒子朝着相对而言质量更大的粒子移动；各粒子的速度和位置从初代开始改变，粒子所在的位置就是问题的可能解，每个粒子的性能与其惯性质量成一定的比例，当算法终止时，选择性能最好的粒子作为问题的解；

所述基于引力搜索算法的MUSIC谱峰搜索方法，将引力搜索算法用于MUSIC算法的谱峰搜索问题，则种群中粒子的位置代表着估计方位角在每次迭代过程中的解，所述方法包括如下步骤：

步骤一，初始化所有粒子在给定维度的速度和位置信息，并设置迭代次数和引力搜索算法中的初始参数，将MUSIC空间谱函数P_MUSIC设为引力搜索算法的适应度函数，将粒子的位置作为自变量θ，代入MUSIC空间谱函数P_MUSIC算出所有粒子的适应度值，并根据适应度值衡量粒子的性能；

步骤二，运用引力搜索算法在每一次迭代过程中更新重力常数，并利用适应度值算出每个粒子的惯性质量，根据重力常数和惯性质量计算出每个粒子的引力之和，再由引力之和和惯性质量得出每个粒子的加速度；

步骤三，利用加速度计算出每次迭代过程中所有粒子的速度和位置，适应度值最大的粒子的位置代表着DOA估计问题的解；

步骤四，检查算法的终止条件，如果满足，终止算法并选择适应度值最大的一组解作为DOA估计问题的解。

2.如权利要求1所述的基于引力搜索算法的MUSIC谱峰搜索方法，其特征在于，所述方法的具体实现过程包括：

步骤1：随机地初始化种群中的粒子，每个粒子都有自己的初始位置，各粒子的位置X_i的表达式为：

其中，

和D分别表示粒子i在第k维搜索空间的位置和搜索空间的维数，D也等于待估计信源的个数，N为初始种群大小，算法在各搜索空间的搜索范围在[-90°,90°]之间；

步骤2：算出各粒子在不同时刻的惯性质量，计算公式如下：

M_ai＝M_pi＝M_i,i＝1,2,..,N (9)

M_ai和M_pi分别代表第i个粒子主动惯性质量和被动惯性质量，M_i(t)代表第i个粒子的惯性质量，t为时刻，fit_i(t)代表t时刻第i个粒子的适应度值，是根据MUSIC空间谱函数

得出的值，其中，a(θ)为信号子空间的方向向量，V_N为噪声子空间，[·]^H指共轭矩阵；

依据上述MUSIC空间谱函数P_MUSIC公式计算每次迭代中所有粒子的适应度，由于该问题是优化最大值问题，时刻t的最优解best(t)和最差解worst(t)分别为：

步骤3：从引力定律中可以推出，第k维上一个粒子i在t时刻受到另一个粒子j的作用力的计算公式如下：

R_ij(t)表示在t时刻两个粒子i和j的欧几里得度量，ε表示一个非常小的常数；

在t时刻的万有引力常数G(t)为：

G(t)＝G₀×e^-α·t/T (18)

G₀表示在t₀时刻G的取值，取值为100，α为下降系数，取值为20，t对应着当前循环次数，T代表最大循环次数；

则在t时刻，X_i在第k维上受到的引力总和是由种群中其余粒子对其产生的全部引力相加得到，计算公式如下所示：

式中，rand_j表示在[0,1]之间的随机数；

步骤4：从牛顿的力学原理可以得知，粒子产生的加速度等于该粒子受到的引力总和与其惯性质量之比，则粒子i的加速度的计算公式如下：

然后在每一次更新过程中，根据计算得到的加速度来更新每个粒子的速度，再根据粒子的速度更新每个粒子的位置，计算公式如下所示：

Images