CN111398890A

CN111398890A - 一种基于协方差矩阵重构的布谷鸟搜索doa估计方法

Info

Publication number: CN111398890A
Application number: CN202010211899.7A
Authority: CN
Inventors: 李鹏; 张伟; 夏翔; 马味敏
Original assignee: Nanjing University of Information Science and Technology
Current assignee: Nanjing University of Information Science and Technology
Priority date: 2020-03-24
Filing date: 2020-03-24
Publication date: 2020-07-10
Anticipated expiration: 2040-03-24
Also published as: CN111398890B

Abstract

本发明公开了一种基于协方差矩阵重构的布谷鸟搜索DOA估计方法，具体包括如下步骤：S1：通过信号接收模型，建立低秩矩阵重构模型；S2：初始化布谷鸟群体数量，随机生成一组鸟窝位置，建立子种群；S3：根据所述低秩矩阵重构模型，确定目标函数，计算每个所述子种群的平均适应度，确定出每个所述子种群的规模；S4：每个所述子种群根据位置更新公式进行鸟窝位置更新；S5：对所述鸟窝位置更新后的子种群进行判定，从中确定出目标信号的DOA估计值。本发明在与对相邻目标信号的分辨能力上，具有更好的均方根误差，并且对于信号环境有更强的适应能力，也能在保证估计精度的前提下，减少计算量，从而达到提升优化速度的效果。

Description

一种基于协方差矩阵重构的布谷鸟搜索DOA估计方法

技术领域

本发明涉及阵列信号处理技术领域，尤其涉及一种基于协方差矩阵重构的布谷鸟搜索DOA估计方法。

背景技术

几十年来，阵列信号处理一直是一个非常活跃的研究领域。作为阵列信号处理中的一个重要研究方向，目标波达方向估计已经广泛应用于水下导航、战场侦察、海上作业等多个方面，引起人们的高度重视。阵列DOA估计也称空间谱估计，其主要研究内容是设计空间阵列结构模型，根据阵列接收到的观测数据，选择高效的、性能优良的算法处理空间内信号的波达角的估计问题。

阵列DOA估计方法早期通过波束形成来实现，主要包括常规波束形成方法和最小方差无畸变响应波束形成方法，这类方法实现简单，应用广泛，但对多个目标信号或者角度相邻的目标信号的分辨能力不佳。现代高精度DOA估计方法主要包括子空间类方法和稀疏重构方法。子空间类方法包括多重信号分类方法和旋转不变子空间方法等；稀疏重构方法包括L1-SVD和稀疏贝叶斯学习等。

子空间类方法可实现对目标信号的高分辨率估计，但其在小快拍、低信噪比条件下的估计性能出现严重恶化的情况，对于信号环境的适应能力较差。在基于稀疏重构的DOA估计方法中，无格稀疏DOA估计方法逐渐成为研究热点。这类方法最大的特点是无需将空间角度离散为网格而直接求解DOA估计值，因此可以避免因为网格划分与真实DOA不匹配所造成的估计误差。文献"A fast gridless covariance matrix reconstruction method forone-and two-dimensional direction-of-arrival estimation"提出了一种无网格的二维测向方法即协方差矩阵重构法。该方法通过低秩矩阵重构模型实现阵列协方差矩阵的恢复，并利用范德蒙分解定理得到DOA的估计值，从根本上避免了网格划分，实现对目标信号的方位角和俯仰角自动配对。由于对于以矩阵秩为目标函数的半定规划问题的最优解是通过CVX工具箱来求解的，算法仍存在计算准确性不足、效率不高的问题。

发明内容

发明目的：针对现有技术存在的对于目标信号DOA估计准确性不高、相邻信号分辨能力不佳、适应信号环境能力较差的问题，本发明提出一种基于协方差矩阵重构的布谷鸟搜索DOA估计方法。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：

一种基于协方差矩阵重构的布谷鸟搜索DOA估计方法，所述布谷鸟搜索DOA估计方法具体包括如下步骤：

S1：通过信号接收模型，建立低秩矩阵重构模型；

S2：初始化布谷鸟群体数量，随机生成一组鸟窝位置，建立子种群；

S3：根据所述低秩矩阵重构模型，确定目标函数，计算每个所述子种群的平均适应度，确定出每个所述子种群的规模；

S4：每个所述子种群根据位置更新公式进行鸟窝位置更新；

S5：对所述鸟窝位置更新后的子种群进行判定，从中确定出目标信号的DOA估计值。

进一步地讲，在所述步骤S1中，建立所述低秩矩阵重构模型，具体如下：

S1.1：通过所述信号接收模型，确定入射信号的俯仰角、入射信号的方位角、入射信号和信号接收模型X轴之间的夹角、入射信号和信号接收模型Y轴之间的夹角四者之间的关系，具体为：

其中：α_k为第k个入射信号和信号接收模型X轴之间的夹角，β_k为第k个入射信号和信号接收模型Y轴之间的夹角，φ_k为第k个入射信号的俯仰角，θ_k为第k个入射信号的方位角；

S1.2：根据所述入射信号的俯仰角、入射信号的方位角、入射信号和信号接收模型X轴之间的夹角、入射信号和信号接收模型Y轴之间的夹角四者之间的关系，获取信号接收模型的输出，具体为：

X＝AS+V

其中：X为信号接收模型的输出，A为信号接收模型的流形矩阵，S为入射目标信号，V为观测噪声；

S1.3：根据所述信号接收模型的输出，得到信号接收模型接收信号的协方差矩阵，具体为：

R＝E[XX^H]

＝AR_sA^H+σI

＝T(u,v)+σI

其中：R为信号接收模型接收信号的协方差矩阵，X为信号接收模型的输出，R_s为目标信号的协方差矩阵，A为信号接收模型的流形矩阵，σ为噪声功率，I为单位矩阵，T(u,v)为阵列的无噪协方差矩阵；

S1.4：对所述阵列的无噪协方差矩阵进行多维范德蒙德分解，同时所述分解后的阵列的无噪协方差矩阵，利用二维CMRA模型得到低秩矩阵重构模型，具体为：

T(u,v)≥0

其中：

T(u,v)为阵列的无噪协方差矩阵，

为样本协方差矩阵，σ为噪声功率，N为信号接收模型中阵元的总个数。

进一步地讲，在所述步骤S2中，建立所述子种群，具体如下：

S2.1：初始化布谷鸟群体数量，根据目标信号的数目，确定子种群的数目，具体为：

M＝K+1

其中：M为子种群的数目，K为目标信号的数目；

S2.2：对所述子种群的生存规模进行最大限制和最小限制；

S2.3：根据所述子种群的数目、子种群生存规模的最大限制和最小限制，确定初始鸟窝位置分布，具体为：

E_i＝E_min+rand(0,1)(E_max-E_min)

其中：E_i为第i个包含α_k和β_k角度信息的二维参数解，E_min为搜索范围的下限值，E_max为搜索范围的上限值，rand(0,1)为0到1之间服从均匀分布的一个随机数，α_k为第k个入射信号和信号接收模型X轴之间的夹角，β_k为第k个入射信号和信号接收模型Y轴之间的夹角。

进一步地讲，所述子种群生存规模的最大限制和最小限制，具体为：

所述子种群生存规模的最大限制为不能超过种群规模的一半；

所述子种群生存规模的最小限制为种群规模的1/10。

进一步地讲，在所述步骤S3中，确定出每个所述子种群的规模，具体如下：

S3.1：根据所述低秩矩阵重构模型，确定目标函数，并将所述低秩矩阵重构模型的函数值作为目标函数的适应度值，所述目标函数具体为：

f(α,β)＝tr[T(u,v)]

其中：f(α,β)为目标函数，T(u,v)为阵列的无噪协方差矩阵；

S3.2：根据所述目标函数，计算出所述种群中每个个体的适应度值，并将所述每个个体的适应度值进行排序；

S3.3：在每个所述子种群中，确定出所有个体的最小适应度值，所述最小适应度值即为子种群的最优解，所述最小适应度值对应的个体位置即为子种群的最优位置；

在所述种群中，确定出所有个体的最小适应度值，所述最小适应度值即为种群的最优解，所述最小适应度值对应的个体位置即为种群的最优位置；

S3.4：计算所述每个子种群的平均适应度，确定出子种群规模，所述平均适应度和子种群规模的求取公式，具体为：

其中：

为第k个子种群第t代的平均适应度值，fit_ki(t)为第k个子种群第t代中第i个个体的适应度值，n_k(t)为第k个子种群第t代的规模，n_k(t+1)为第k个子种群第t+1代的规模，K为目标信号的数目，Q为初始化布谷鸟群体数量。

进一步地讲，在所述步骤S4中，每个所述子种群根据位置更新公式进行鸟窝位置更新，所述更新后的位置具体为：

其中：E_k(t+1)为第k个子种群第t+1次迭代的位置，E_k(t)为第k个子种群第t次迭代的位置，ρ为惯性权重，E_kbest(t)为第k个子种群第t次迭代的最优位置，Levy(ω)为莱维飞行模式。

进一步地讲，在所述步骤S5中，确定出所述目标信号的DOA估计值，具体如下：

S5.1：根据所述鸟窝位置更新后的子种群对应的最优解和最优位置，将所述鸟窝位置更新后的子种群进行同种互斥判定；

S5.2：对在进行所述同种互斥判定过程中，已经进化停滞的子种群进行新老更替，并产生一组等同规模的新解，同时获取所述新产生的子种群中的每个个体的自适应度值，并继续进行进化；

S5.3：当迭代次数达到预设的最大迭代次数或算法满足收敛条件时，停止搜索，并结束进化过程，若未满足，则返回步骤S4，当所有的子种群均满足收敛性判定时，则输出最优解和最优位置；

S5.4：所述输出的最优解和最优位置根据范德蒙分解定理，确定出所述目标信号的DOA估计值。

进一步地讲，在所述步骤S5.1中，将所述鸟窝位置更新后的子种群进行同种互斥判定，具体如下：

S5.1.1：获取两个不同子种群之间的准则值，所述准则值的计算公式，具体为：

其中：I_ij为第i个子种群和第j个子种群之间的准则值，fit_i(t)为第i个子种群的最优解适应度值，fit_j(t)为第j个子种群的最优解适应度值，

为第i个子种群和第j个子种群之间最优解位置参数的欧氏距离；

S5.1.2：将所述两个不同子种群之间的准则值和预先设定的阈值进行比较，当所述准则值小于预先设定的阈值，则所述两个不同子种群搜索到相似的最优解位置，反之，则所述两个不同子种群未搜索到相似的最优解位置；

S1.5.3：将所述搜索到相似的最优解位置的两个不同子种群的平均适应度值进行比较，从中选出较大平均适应度值对应的子种群，并将所述较大平均适应度值对应的子种群进行删除。

进一步地讲，在对所述已经进化停滞的子种群进行新老更替的过程中，其中包含有最优个体的子种群需要进行保留。

进一步地讲，在对已经进化停滞的子种群进行新老更替，并产生一组等同规模的新解的过程中，所述新解的确定公式和初始鸟窝位置分布的确定公式相同。

有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：

本发明相较于2-D CMRA算法，在与对相邻目标信号的分辨能力上，具有更好的均方根误差，并且对于信号环境有更强的适应能力，也能在保证估计精度的前提下，减少计算量，从而达到提升优化速度的效果，同时本发明相较于其他优化方法，具有估计精度高、优化速度快、兼顾局部与全局搜索能力等特点。

附图说明

图1是本发明布谷鸟搜索DOA估计方法的流程示意图；

图2是本发明的阵列模型图；

图3是入射信号角度间隔较大时的DOA估计散点图；

图4是入射信号角度间隔较小时的DOA估计散点图；

图5是DOA估计RMSE随信噪比SNR的变化对比图；

图6是DOA估计RMSE随快拍数的变化对比图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。其中，所描述的实施例是本发明一部分实施例，而不是全部的实施例。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。

实施例1

参考图1，本实施例提供了一种基于协方差矩阵重构的布谷鸟搜索DOA估计方法，该布谷鸟搜索DOA估计方法首先基于协方差匹配准则，通过低秩矩阵重构模型实现阵列协方差矩阵的恢复，由范德蒙分解定理建立一个迹最小化问题。其次改进布谷鸟搜索路径的多个参数，并将隔离小生境技术应用于布谷鸟搜索算法，求解出迹最小化问题在空间内的多个最优解，从而对多个入射目标的二维角度进行估计。具体过程如下：

步骤S1：通过信号接收模型，建立低秩矩阵重构模型。本实施例主要针对均匀矩形面阵进行说明，如图2所示。在整个面阵中一共有N个阵元，其中沿面阵的X轴上的阵元数为N_x个，沿面阵的Y轴上的阵元数为N_y个，各阵元之间的间隔为d，入射信号源的波长为λ。同时沿面阵的X轴上的阵元数N_x、沿面阵的Y轴上的阵元数N_y和整个面阵总阵元数N之间满足如下关系，具体为：

N＝N_x×N_y

其中：N为整个面阵总阵元数，N_x为沿面阵的X轴上的阵元数，N_y为沿面阵的Y轴上的阵元数。

各阵元之间的间隔d和入射信号源的波长λ之间满足如下关系，具体为：

d＝λ2

其中：d为各阵元之间的间隔，λ为入射信号源的波长。

在本实施例中，建立低秩矩阵重构模型，具体如下：

步骤S1.1：本实施例中，有K个远场窄带信号入射到该均匀矩形面阵上，其中第k个入射信号的俯仰角为φ_k和方位角为θ_k，第k个入射信号和该均匀矩形面阵X轴之间的夹角为α_k，第k个入射信号和该均匀矩形面阵Y轴之间的夹角为β_k，同时φ_k、θ_k、α_k和β_k四者之间的关系，具体为：

其中：α_k为第k个入射信号和信号接收模型X轴之间的夹角，β_k为第k个入射信号和信号接收模型Y轴之间的夹角，φ_k为第k个入射信号的俯仰角，θ_k为第k个入射信号的方位角。

从而可以得到如下关系式，具体为：

从上式可以看出俯仰角φ_k、方位角θ_k与和信号接收模型X轴之间的夹角α_k、信号接收模型Y轴之间的夹角β_k之间存在有一定的转换关系，因此可以通过信号接收模型X轴之间的夹角α_k、信号接收模型Y轴之间的夹角β_k确定出入射信号的角度信息。

步骤S1.2：根据步骤S1中入射信号的俯仰角φ_k、入射信号的方位角θ_k、入射信号和信号接收模型X轴之间的夹角α_k、入射信号和信号接收模型Y轴之间的夹角β_k四者之间的关系，当阵列接收到L个快拍时，可以确定出信号接收模型的输出，具体为：

X＝AS+V

其中：X为信号接收模型的输出，A为信号接收模型的流形矩阵，S为入射目标信号，V为观测噪声。

具体地讲，信号接收模型的输出X，具体为：

X＝[x(1),…,x(L)]

其中：X为信号接收模型的输出，x(1)为第一时刻阵列输出的矢量化形式，x(L)为第L时刻阵列输出的矢量化形式。

信号接收模型的流形矩阵A，具体为：

A＝A_x⊙A_y

其中：A为信号接收模型的流形矩阵，A_x为信号接收模型X轴上的均匀线性阵列的阵列流形矩阵，A_y为信号接收模型Y轴上的均匀线性阵列的阵列流形矩阵。

同时信号接收模型X轴上的均匀线性阵列的阵列流形矩阵A_x、信号接收模型Y轴上的均匀线性阵列的阵列流形矩阵A_y，具体为：

其中：

A_x为信号接收模型X轴上的均匀线性阵列的阵列流形矩阵，A_y为信号接收模型Y轴上的均匀线性阵列的阵列流形矩阵，a_x(α_K)为信号接收模型X轴上的均匀线性阵列关于来向α_K的导向矢量，α_k为第k个入射信号和信号接收模型X轴之间的夹角，a_y(β_K)为信号接收模型Y轴上的均匀线性阵列关于来向β_K的导向矢量，β_k为第k个入射信号和信号接收模型Y轴之间的夹角，N_x为沿面阵的X轴上的阵元数，N_y为沿面阵的Y轴上的阵元数。

入射目标信号S，具体为：

S＝[s(1),…,s(L)]

其中：S为入射目标信号，s(1)为第一时刻入射目标信号的矢量化形式，s(L)为第L时刻入射目标信号的矢量化形式。

观测噪声V，具体为：

V＝[v(1),…,v(L)]

其中：V为观测噪声，v(1)为第一时刻观测噪声的矢量化形式，v(L)为第L时刻观测噪声的矢量化形式。

步骤S1.3：根据步骤S1.2中得到的信号接收模型的输出，可以得到信号接收模型接收信号的协方差矩阵，具体为：

R＝E[XX^H]

＝AR_sA^H+σI

＝T(u,v)+σI

其中：R为信号接收模型接收信号的协方差矩阵，X为信号接收模型的输出，R_s为目标信号的协方差矩阵，A为信号接收模型的流形矩阵，σ为噪声功率，I为单位矩阵，T(u,v)为阵列的无噪协方差矩阵。

具体地讲，阵列的无噪协方差矩阵T(u,v)是一个嵌套Toeplitz矩阵，可以表示如下形式，具体为：

其中：T(u,v)为阵列的无噪协方差矩阵，A为信号接收模型的流形矩阵，S为入射目标信号，a_x(α_K)为信号接收模型X轴上的均匀线性阵列关于来向α_K的导向矢量，α_k为第k个入射信号和信号接收模型X轴之间的夹角，a_y(β_K)为信号接收模型Y轴上的均匀线性阵列关于来向β_K的导向矢量，β_k为第k个入射信号和信号接收模型Y轴之间的夹角，p_k为阈值。

从上式可以看出，阵列的无噪协方差矩阵T(u,v)是一个半正定矩阵，具有低秩性质，且具有入射信号和信号接收模型X轴之间的夹角α_k、入射信号和信号接收模型Y轴之间的夹角β_k的角度信息。从而可以通过对阵列的无噪协方差矩阵T(u,v)进行多维范德蒙德分解，得到二维角度的估计值。

具体地讲，多维范德蒙德分解定理具体为：

给定一个秩为r＜min{N_x,N_y}的半正定嵌套Toeplitz矩阵T(u,v)，使其满足如下分解：

同时p_k＞0，且该分解具有唯一性。

其中：T(u,v)为阵列的无噪协方差矩阵，a_x(α_K)为信号接收模型X轴上的均匀线性阵列关于来向α_K的导向矢量，α_k为第k个入射信号和信号接收模型X轴之间的夹角，a_y(β_K)为信号接收模型Y轴上的均匀线性阵列关于来向β_K的导向矢量，β_k为第k个入射信号和信号接收模型Y轴之间的夹角，p_k为阈值。

步骤S1.4：为了进一步求解出阵列的无噪协方差矩阵T(u,v)，利用二维CMRA模型得到低秩矩阵重构模型，具体为：

T(u,v)≥0

其中：

T(u,v)为阵列的无噪协方差矩阵，

具体地讲，样本协方差矩阵

在实际阵列中采样点数是有限的，故取协方差矩阵估计式。具体为：

其中：

为样本协方差矩阵，X为信号接收模型的输出，L为信号接收模型接收到快拍个数。

当得到低秩矩阵重构模型最小化问题的最优解后，即可利用多维范德蒙德分解定理对入射角度进行估计。

步骤S2：现有技术对于以矩阵迹为目标函数的半定规划问题的最优解是通过CVX工具箱来求解的，仍存在计算准确性不足、效率不高的问题。因此本实施例中采用一种基于隔离小生境的布谷鸟搜索算法对问题进行优化后求解出需要的多个最优解。

根据目标信号的数目K，将初始化布谷鸟群体数量设置为Q，随机生成一组鸟窝位置，建立子种群，具体如下：

步骤S2.1：将初始化布谷鸟群体数量设置为Q，根据目标信号的数目K，确定出子种群的数目M，目标信号的数目K和子种群的数目M之间满足如下关系，具体为：

M＝K+1

其中：M为子种群的数目，K为目标信号的数目。

步骤S2.2：对子种群的生存规模进行最大限制和最小限制，具体为：

子种群生存规模的最大限制为不能超过种群规模的一半，即

M_max＝Q/2

其中：M_max为子种群的最大数目，Q为初始化布谷鸟群体数量。

子种群生存规模的最小限制为种群规模的1/10，即

M_min＝Q/10

其中：M_min为子种群的最小数目，Q为初始化布谷鸟群体数量。

步骤S2.3：根据子种群的数目M、子种群生存规模的最大限制M_max和最小限制M_min，可以确定出初始鸟窝位置分布，具体为：

E_i＝E_min+rand(0,1)(E_max-E_min)

步骤S3：根据步骤S1.4中得到的低秩矩阵重构模型，确定出目标函数，并计算每个子种群的平均适应度，确定出每个子种群的规模。具体如下：

步骤S3.1：根据步骤S1.4中得到的低秩矩阵重构模型，确定出目标函数，并将低秩矩阵重构模型的函数值作为目标函数的适应度值，该目标函数具体为：

f(α,β)＝tr[T(u,v)]

其中：f(α,β)为目标函数，T(u,v)为阵列的无噪协方差矩阵。

步骤S3.2：根据步骤S3.1中的目标函数f(α,β)，计算出种群中每个个体的适应度值，并将每个个体的适应度值进行排序，排序可以按照从大到小的顺序进行排列，也可以按照从小到大的顺序进行排列。

步骤S3.3：在每个子种群中，将该种群中所有个体的应度值进行比较，从中确定出该种群中所有个体的最小适应度值，具体地讲，该最小适应度值即为子种群的最优解，该最小适应度值对应的个体位置即为子种群的最优位置。

在整个种群中，将整个种群中所有个体的应度值进行比较，从中确定出整个种群中所有个体的最小适应度值，具体地讲，该最小适应度值即为整个种群的最优解，该最小适应度值对应的个体位置即为整个种群的最优位置。

步骤S3.4：根据步骤S3.2中得到的每个个体的适应度值，计算出每个子种群的平均适应度，并确定出子种群规模，该平均适应度和子种群规模的求取公式，具体为：

其中：

步骤S4：各个子种群均根据位置更新公式进行鸟窝位置更新，该更新后的位置具体为：

在本实施例中，惯性权重ρ具体为：

ρ＝0.158-0.058exp((t/t_max)^1.5)

其中：ρ为惯性权重，t为当前的迭代次数，t_max为最大迭代次数。

莱维飞行模式Levy(ω)具体为：

w和z均服从正态分布，具体为：

ω采用自适应策略，具体为：

ω＝0.582exp((t/t_max)^1.8)+0.218

其中：

Levy(ω)为莱维飞行模式，t为当前的迭代次数，t_max为最大迭代次数。

从上式可以看出，在算法前期，ω的取值相对较小，相应的莱维飞行模式Levy(ω)值较大，因此全局搜索能力较强。在在算法后期，ω的取值相对较大，相应的莱维飞行模式Levy(ω)值较小，因此局部搜索能力较强

这里的惯性权重ρ和ω的自适应取值策略，能让算法保持很强的全局搜索寻优能力，同时兼顾一定的局部精细搜索能力。

步骤S5：对鸟窝位置更新后的子种群进行判定，从中确定出目标信号的DOA估计值，具体如下：

步骤S5.1：根据鸟窝位置更新后的子种群对应的最优解和最优位置，将鸟窝位置更新后的子种群进行同种互斥判定，即当两个不同子种群搜索到相似的最优解位置，则将平均适应度值较大的子种群舍弃掉，并重新生成一组同等规模的新解取代旧的解。

具体如下：

步骤S5.1.1：获取两个不同子种群之间的准则值，该准则值的计算公式，具体为：

为第i个子种群和第j个子种群之间最优解位置参数的欧氏距离。

步骤S5.1.2：将两个不同子种群之间的准则值和预先设定的阈值进行比较，当得到的准则值小于预先设定的阈值时，则两个不同子种群搜索到相似的最优解位置，反之，则两个不同子种群未搜索到相似的最优解位置。

步骤S1.5.3：将搜索到相似的最优解位置的两个不同子种群的平均适应度值进行比较，从中选出较大平均适应度值对应的子种群，并将较大平均适应度值对应的子种群进行删除。

步骤S5.2：对在进行上述同种互斥判定过程中，已经进化停滞的子种群进行新老更替，即将需要删除的子种群进行新老更替，同时产生一组等同规模的新解，新解的求取和步骤S2.3中初始鸟窝位置分布的确定公式相同。在本实施例中，还要获取新产生的子种群中的每个个体的自适应度值，并对其继续进行进化。

值得注意的是，在对已经进化停滞的子种群进行新老更替的过程中，其中包含有最优个体的子种群需要进行保留。

步骤S5.3：当迭代次数达到预设的最大迭代次数或算法满足收敛条件时，停止搜索，并结束进化过程，若未满足，则返回步骤S4继续运行。当所有的子种群均满足收敛性判定时，则输出最优解和最优位置。

步骤S5.4：输出的最优解和最优位置根据范德蒙分解定理，可以得到入射信号和信号接收模型X轴之间的夹角α_k、入射信号和信号接收模型Y轴之间的夹角β_k，经过步骤S1.1中的变换公式，即可得到入射信号的俯仰角φ_k、入射信号的方位角θ_k，从而即可确定出目标信号的DOA估计值。

在本实施例中，通过以下仿真进一步进行说明，

仿真条件：

考虑在整个均匀矩形阵列中一共有N＝N_x×N_y个阵元，其中N_x＝N_y＝4，各阵元之间的间隔为d＝λ/2，噪声为高斯白噪声。信源个数K＝3，布谷鸟种群数量选择为100，子种群数量选择为4个，最大迭代次数为100。

在对DOA估计方法的估计准确性进行评估时，以各入射信号的平均测向均方根误差作为指标，平均测向均方根误差定义为：

其中：RMSE为平均测向均方根误差，F为实验次数，K为信源个数，

和

为第n次仿真实验中的二维角度估计值，α⁽ⁿ⁾和β⁽ⁿ⁾为第n次仿真实验中的二维角度真实值。

仿真内容：

仿真实验1：考虑三个入射角度分别为α＝[-45°，-20°，5°]和β＝[5°,30°,0°]的窄带非相干信号，快拍数为200，信噪比为20dB，进行50次蒙特卡洛实验，仿真结果如图3所示。由图3可知，从散点集中程度上看，2-D CMRA方法的DOA估计值散点图具有较大的离散程度，而本实施例所提算法的散点图较为集中，具有更加稳定的估计性能。从散点偏离程度上看，本实施例所提算法的估计值散点与真实信号方向的偏移量可以忽略不计，而2-D CMRA方法的色块与真实信号方向的偏移量明显大于本实施例所提算法。由此表明，相较于2-D CMRA方法，本实施例所提方法具有更高的估计精度和更加稳定的估计性能。

仿真实验2：考虑三个入射角度分别为α＝[-20°，-15°，-5°]和β＝[20°,5°,15°]的窄带非相干信号，快拍数为200，信噪比为20dB，进行50次蒙特卡洛实验，仿真结果如图4所示。由图4可知，从散点集中程度上看，2-D CMRA方法的DOA估计值散点图的离散程度进一步扩大，而本实施例所提算法的散点图与图3相比虽略有扩散但仍具有较高的集中程度。从散点偏离程度上看，2-D CMRA算法的色块明显偏离信号真实方向，与图3相比其偏离程度大幅加剧，产生这种偏差的主要原因是几个入射信号方向过于靠近。由此表明，当入射信号角度较为接近时，相较于2-D CMRA方法，本实施例所提算法依然具有更好的稳定性能和更高的估计准确性，同时对相邻入射信号有较强的分辨能力。

仿真实验3：考虑入射信号角度同实验1，快拍数为200，进行50次蒙特卡洛实验，仿真结果如图5所示。由图5可知，当信噪比为-15dB时，本实施例所提算法的RMSE在10以下，低于其他两种算法，表明所提算法在低信噪比条件下具有较好的估计性能；随着信噪比的增加，三种算法的RMSE逐渐降低，估计性能逐渐增强；当信噪比大于10dB时，本实施例所提算法的RMSE小于0.1，且能够较好地逼近CRLB下界，表明相较于其他两种算法，所提算法具有较好的角度估计准确性和适应信号环境能力。

仿真实验4：考虑入射信号角度同实验1，信噪比为10dB，进行50次蒙特卡洛实验，仿真结果如图6所示。由图6可知，随着快拍数的增加，本实施例所提算法、2-D ESPRIT、2-DCMRA算法的DOA估计RMSE均逐渐降低。并且，所提算法的RMSE显著低于其他算法，特别在快拍数较少的条件下更加明显。由此表明，在信噪比相同情况下，本实施例所提算法的DOA估计性能优于其他两种算法，具有较高的准确性和适应信号环境能力。

以上示意性的对本发明及其实施方式进行了描述，该描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际的结构和方法并不局限于此。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均属于本发明的保护范围。