CN113589255A - 一种基于多频联合稀疏贝叶斯学习的到达角估计方法 - Google Patents

Abstract

本发明提出了一种基于多频联合稀疏贝叶斯学习的到达角估计方法。本发明获取多个频率通道的快拍数据；构建频率通道快拍数据对应的稀疏模型，构建利用所有频率通道的快拍数据联合估计到达角的稀疏模型；将频率通道的快拍数据对应的稀疏模型通过多频联合稀疏贝叶斯学习算法估计的目标到达角。本发明优点在于，与已有的多频到达角估计算法相比，该算法具有更好的估计精度和分辨率，尤其是在低信噪比和小快拍场景下。

Description

一种基于多频联合稀疏贝叶斯学习的到达角估计方法

技术领域

本发明属于雷达信号处理领域，具体涉及一种基于多频联合稀疏贝叶斯学习的到达角估计方法。

背景技术

近年来，随着低空领域的逐步开放和低空监测需求的日趋增加，外源雷达受到越来越多的关注。外源雷达常利用双基距离和DOA定位和跟踪目标。因此，精准度高和稳定性强的外源雷达DOA估计方法成为了研究的重点。

针对多频外源雷达的特点，本发明提出了一种基于MFJSBL的DOA估计算法。该算法利用目标在空域的稀疏性以提高估计精度和分辨率。首先对提取的频率通道数据通过杂波抑制、相干积累、波束形成、恒虚警率检测等一系列预处理，获取快拍数据。然后，在稀疏贝叶斯学习框架下，联合提取多个频率通道的数据，通过使证据函数最大来估计目标DOA。仿真结果表明，与现有的多频DOA估计算法相比，该算法具有更好的估计精度和分辨率，尤其是在低信噪比和小快拍场景下。最后，本发明采用基于多个调频广播的外源雷达采集的外场实测数据验证了提出方法的有效性。

发明内容

针对多频外辐射源雷达的特点，本发明提出了一种基于多频联合稀疏贝叶斯学习的到达角估计方法。

本发明的技术方案为一种基于多频联合稀疏贝叶斯学习的到达角估计方法，具体包含以下步骤：

步骤1：获取多个频率通道的快拍数据；

步骤2：构建步骤1所述的第n个频率通道快拍数据对应的稀疏模型，构建利用所有频率通道的快拍数据联合估计到达角的稀疏模型；

步骤3：将第n个频率通道的快拍数据对应的稀疏模型通过多频联合稀疏贝叶斯学习算法估计的目标到达角；

作为优选，所述步骤1具体实现包括以下子步骤：

步骤1.1：构建第t个时刻采集的第m个监测通道接收到的第n个频率通道的信号，具体如下：

式中，

表示第t个时刻采集的第n个频率通道的参考信号；

为直达波信号相对于参考信号的归一化幅度；

为直达波相对于参考信号的时延；

为直达波的到达角；

表示第p条多径信号相对于参考信号的复幅度；

表示第p条多径信号相对于参考信号的时延；

表示第p条多径信号的DOA；η_k表示第k个目标的幅度；τ_k表示第k个目标的时延；

分别表示第k个目标的多普勒频率；θ_k表示第k个目标的DOA；

表示第n个频率通道接收的第m个监测信号中的高斯白噪声；

表示监测信号中第n个频率通道的目标、直达波和多径信号的阵列导向矢量；M表示监测通道数，N表示频率通道数，K表示目标数，N_c表示多径数，T_int＝N_intT_s表示每场数据积累时间，N_int表示积累数据点数，T_s＝1/F_s表示离散时间间隔，F_s表示采样率；

步骤1.2：将第t个时刻采集的第m个监测通道接收到的第n个频率通道的信号利用杂波抑制算法滤除监测信号中明显强于目标的直达波和多径杂波，得到抑制后第t个时刻第m个监测通道接收的第n个频率通道的信号为：

步骤1.3：通过第m个监测通道接收的第n个频率通道的信号与相应的参考信号互相关生成的距离多普勒谱为：

τ＝[0,T_s,...,(R_max/ΔR)*T_s]；f_d＝[(-v_max/Δv)*Δf,...,0,...,(v_max/Δv)*Δf]

式中，(·)^*为共轭运算符；R_max表示外源雷达可探测的最大距离；v_max表示外源雷达可探测的最大速度；ΔR＝c/B表示可分辨的距离间隔，c表示光速，B表示照射源信号带宽；Δv＝λ*Δf表示可分辨的速度间隔，λ表示照射源信号波长，Δf＝1/T_int表示可分辨的多普勒间隔；

步骤1.4：利用恒虚警率检测方法获取第m个监测通道接收的第n个频率通道的信号与相应的参考信号互相关生成的距离多普勒谱中目标的复幅度值作为到达角估计的快拍数据；假定K个目标位于同一距离多普勒，则第m个监测通道接收的第n个频率通道的信号与相应的参考信号互相关生成的距离多普勒谱中第k个时延τ_k和第k个多普勒

处的目标幅度值为：

式中，

表示第n个频率通道接收的第k个目标在时延和多普勒

处的自模糊函数值；

表示第n个频率通道接收的第m个监测信号相干积累后的噪声；

步骤1.5：整理估计目标DOA所需的快拍数据。这里需指出的，步骤1.4直接获取的第n个频率通道的快拍数据为

是单个的快拍数据；

其中，

进一步地，为增加估计精度，采用频域滑窗的方式计算步骤1.4的距离多普勒谱L次，得到L个快拍数据；

整理第n个频率通道的L个快拍数据为矩阵形式为：

Y_n＝A_nC_n+Z_n n＝1,...,N

式中，

Y_n＝[y^n,1,...,y^n,l,...,y^n,L]，

A_n＝[a_n(θ₁),...,a_n(θ_k),...,a_n(θ_K)]，

Y_n表示第n个频率通道对应的所有监测通道的快拍数据；y^n,1表示第n个频率通道的第1快拍在所有监测通道上的数据；y^n,l表示表示第n个频率通道的第l快拍在所有监测通道上的数据；y^n,L表示表示第n个频率通道的第L快拍在所有监测通道上的数据；

表示第n个频率通道的第1个监测通道上获取的第l个快拍数据；

表示第n个频率通道的第m个监测通道上获取的第l个快拍数据；

表示第n个频率通道的第M个监测通道上获取的第l个快拍数据；上标T为转置运算符，L为快拍数。A_n表示第n个频率通道接收的所有目标对应的阵列流型矩阵；a_n(θ₁)表示第n个频率通道接收的第1个目标对应的阵列流型矢量；a_n(θ_k)表示第n个频率通道接收的第k个目标对应的阵列流型矢量；a_n(θ_K)表示第n个频率通道接收的第K个目标对应的阵列流型矢量；

表示第n个频率通道接收的第k个目标在第1个监测通道的阵列流型；

表示第n个频率通道接收的第k个目标在第m个监测通道的阵列流型；

表示第n个频率通道接收的第k个目标在第M个监测通道的阵列流型；C_n为第n个频率通道接收目标的复幅度；Z_n为第n个频率通道的高斯白噪声。

作为优选，步骤2所述构建步骤1所述的第n个频率通道快拍数据对应的基于网格的稀疏模型为：

Y_n＝Φ_nW_n+Z_n n＝1,...,N

式中，

W_n＝[w_n,1,...,w_n,l,...,w_n,L]；Φ_n为第n个频率通道的过完备字典矩阵；

为第1个网格角

对应的第n个频率通道的阵列流型向量；

为第k个网格角

对应的第n个频率通道的阵列流型向量；

为第

个网格角

对应的第n个频率通道的阵列流型向量；W_n为第n个频率通道的过完备字典矩阵对应的所有快拍数据的系数矩阵；w_n,1为第n个频率通道的第1个快拍数据对应的系数向量；w_n,l为第n个频率通道的第l个快拍数据对应的系数向量；w_n,L为第n个频率通道的第L个快拍数据对应的系数向量；

为整个探测角度离散化的网格点数；

作为优选，所述步骤3实现包括以下子步骤：

步骤3.1：构建系数矩阵

服从的概率分布模型为：

式中，

为所有频率通道的系数矩阵；W₁为第1个频率通道对应的系数矩阵；W_n为第n个频率通道对应的系数矩阵；W_N为第N个频率通道对应的系数矩阵；Γ＝diag(γ)为协方差矩阵；diag(γ)表示以γ作为对角元素的对角矩阵；

γ_i≥0为方差，又称为超参数先验，控制着模型的稀疏性；

表示复高斯分布符号。MFJSBL算法假定所有系数矩阵共享超参数等同于为系数矩阵放置相同的稀疏约束，以增加DOA的拟合概率。

步骤3.2：估计超参数先验γ；

在稀疏贝叶斯学习框架下，依据第Ⅱ类最大似然法，通过使证据函数

相对于γ等于零，并引入因子

得到第i个超参数的更新迭代方程为：

式中，

和

表示所有频率通道的快拍数据；Y₁表示第1个频率通道的快拍数据；Y_n表示第n个频率通道的快拍数据；Y_N表示第N个频率通道的快拍数据；β表示所有频率通道的噪声方差；

表示第1个频率通道的噪声方差；

表示第n个频率通道的噪声方差；

表示第N个频率通道的噪声方差；R_n为第n个频率通道的采样协方差矩阵；

为第n个频率通道的协方差矩阵；

为原先的第i个超参数；

为更新的第i个超参数；a_n,i为第n个频率通道的第i个网格对应的方向系数；上标H是共轭转置运算符，tr(·)表示矩阵的迹；

步骤3.3：估计噪声方差β；

同样依据第Ⅱ类最大似然法，通过使证据函数

的导数相对于β的导数等于零得到第n个频率通道的噪声方差

的更新方程为

式中，

为由Φ_n中与γ中η个最大的峰对应的原子向量组成的索引集

构成。如果目标数K是已知的，峰的数目η＝K；否则，估计的目标数依据Akaike信息准则得到。

步骤3.4：计算算法收敛率；

算法收敛率∈定义为超参数γ的相对改变率，写为

当算法收敛率达到预设的值时停止计算。

步骤3.5：依据得到的超参数γ得到估计的目标到达角；

算法实现时，需首先等间隔离散整个探测区域为

的角度网格。超参数γ值离散的角度网格相对应，因此，依据估计的γ值可得到目标的到达角；

本发明优点在于，与已有的多频DOA估计算法相比，该算法具有更好的估计精度和分辨率，尤其是在低信噪比和小快拍场景下。

附图说明

图1：是宽带外源雷达工作原理示意图

图2：是本发明的实现流程图

图3：是均匀线阵(ULA)在快拍数L＝15，探测θ＝0°的单目标，频率通道数为n＝1,3,5时，DOA估计的均方根误差(RMSE)与目标积累SNR的关系。

图4：是ULA在快拍数L＝15，探测θ＝(0°,5°)或θ＝(0°,50°)的两目标，频率通道数为n＝1,3,5时，DOA估计的RMSE与目标积累SNR的关系。

图5：是ULA在快拍数L＝1，探测θ＝0°单目标，频率通道数为n＝1,3,5时，DOA估计的RMSE与目标积累SNR的关系。

图6：是均匀圆形阵列(UCA)在快拍数L＝1，探测θ＝0°单目标，频率通道数为n＝1,3,5时，DOA估计的RMSE与目标积累SNR的关系。

图7：是各种DOA估计方法的RMSE

图8：是DOA估计的RMSE与两目标角度间隔

图9：是单快拍数据估计多目标的归一化功率谱

具体实施方式

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

图1是宽带外源雷达工作原理示意图。图2是本发明的实现流程图。本发明在科学计算软件matlab R2018仿真平台对发明实例予以验证。仿真基于FM广播的外源雷达，其监测天线采用阵元间隔为1.5米的7元均匀线阵。当频率通道数为1，3和5时，其相应的载波频率集合分别为{97.4}M，{97.4,101.8,102.6}MHz，{97.4,98.6,101.8,102.6,103.7}MHz。本发明从目标数、快拍数和阵列类型等方面研究了提出算法对DOA估计精度的影响。

下面结合图1至图9论述本发明的实施，本发明实施的具体步骤为：

步骤1：获取多个频率通道的快拍数据。

具体实现包括以下子步骤：

步骤1.1：构建第t个时刻采集的第m个监测通道接收到的第n个频率通道的信号，具体如下：

式中，

表示第t个时刻采集的第n个频率通道的参考信号；

为直达波信号相对于参考信号的归一化幅度；

为直达波相对于参考信号的时延；

为直达波的到达角；

表示第p条多径信号相对于参考信号的复幅度；

表示第p条多径信号相对于参考信号的时延；

表示第p条多径信号的DOA；η_k表示第k个目标的幅度；τ_k表示第k个目标的时延；

分别表示第k个目标的多普勒频率；θ_k表示第k个目标的DOA；

表示第n个频率通道接收的第m个监测信号中的高斯白噪声；

表示监测信号中第n个频率通道的目标、直达波和多径信号的阵列导向矢量；M＝7表示监测通道数，N＝1,3,5表示频率通道数，K＝1,2表示目标数，N_c＝9表示多径数，T_int＝N_intT_s＝1s表示每场数据积累时间，N_int表示积累数据点数，T_s＝1/F_s表示离散时间间隔，F_s＝300KHz表示采样率；

步骤1.2：将第t个时刻采集的第m个监测通道接收到的第n个频率通道的信号利用杂波抑制算法滤除监测信号中明显强于目标的直达波和多径杂波，得到抑制后第t个时刻第m个监测通道接收的第n个频率通道的信号为：

步骤1.3：通过第m个监测通道接收的第n个频率通道的信号与相应的参考信号互相关生成的距离多普勒谱为：

τ＝[0,T_s,...,(R_max/ΔR)*T_s]；f_d＝[(-v_max/Δv)*Δf,...,0,...,(v_max/Δv)*Δf]

式中，(·)^*为共轭运算符；R_max表示外源雷达可探测的最大距离；v_max表示外源雷达可探测的最大速度；ΔR＝c/B表示可分辨的距离间隔，c表示光速，B表示照射源信号带宽；Δv＝λ*Δf表示可分辨的速度间隔，λ表示照射源信号波长，Δf＝1/T_int表示可分辨的多普勒间隔；

步骤1.4：利用恒虚警率检测方法获取第m个监测通道接收的第n个频率通道的信号与相应的参考信号互相关生成的距离多普勒谱中目标的复幅度值作为到达角估计的快拍数据；假定K个目标位于同一距离多普勒，则第m个监测通道接收的第n个频率通道的信号与相应的参考信号互相关生成的距离多普勒谱中第k个时延τ_k和第k个多普勒

处的目标幅度值为：

式中，

表示第n个频率通道接收的第k个目标在时延和多普勒

处的自模糊函数值；

表示第n个频率通道接收的第m个监测信号相干积累后的噪声；

步骤1.5：整理估计目标DOA所需的快拍数据。这里需指出的，步骤1.4直接获取的第n个频率通道的快拍数据为

是单个的快拍数据；

其中，

进一步地，为增加估计精度，采用频域滑窗的方式计算步骤1.4的距离多普勒谱L次，得到L个快拍数据；

整理第n个频率通道的L个快拍数据为矩阵形式为：

Y_n＝A_nC_n+Z_n n＝1,...,N

式中，

Y_n＝[y^n,1,...,y^n,l,...,y^n,L]，

A_n＝[a_n(θ₁),...,a_n(θ_k),...,a_n(θ_K)]，

Y_n表示第n个频率通道对应的所有监测通道的快拍数据；y^n,1表示第n个频率通道的第1快拍在所有监测通道上的数据；y^n,l表示表示第n个频率通道的第l快拍在所有监测通道上的数据；y^n,L表示表示第n个频率通道的第L快拍在所有监测通道上的数据；

表示第n个频率通道的第1个监测通道上获取的第l个快拍数据；

表示第n个频率通道的第m个监测通道上获取的第l个快拍数据；

表示第n个频率通道的第M个监测通道上获取的第l个快拍数据；上标T为转置运算符，L为快拍数。A_n表示第n个频率通道接收的所有目标对应的阵列流型矩阵；a_n(θ₁)表示第n个频率通道接收的第1个目标对应的阵列流型矢量；a_n(θ_k)表示第n个频率通道接收的第k个目标对应的阵列流型矢量；a_n(θ_K)表示第n个频率通道接收的第K个目标对应的阵列流型矢量；

表示第n个频率通道接收的第k个目标在第1个监测通道的阵列流型；

表示第n个频率通道接收的第k个目标在第m个监测通道的阵列流型；

表示第n个频率通道接收的第k个目标在第M个监测通道的阵列流型；C_n为第n个频率通道接收目标的复幅度；Z_n为第n个频率通道的高斯白噪声。

步骤2：构建步骤1所述的第n个频率通道快拍数据对应的稀疏模型，构建利用所有频率通道的快拍数据联合估计到达角的稀疏模型。

步骤2所述构建步骤1所述的第n个频率通道快拍数据对应的基于网格的稀疏模型为：

Y_n＝Φ_nW_n+Z_n n＝1,...,N

式中，

W_n＝[w_n,1,...,w_n,l,...,w_n,L]；Φ_n为第n个频率通道的过完备字典矩阵；

为第1个网格角

对应的第n个频率通道的阵列流型向量；

为第k个网格角

对应的第n个频率通道的阵列流型向量；

为第

个网格角

对应的第n个频率通道的阵列流型向量；W_n为第n个频率通道的过完备字典矩阵对应的所有快拍数据的系数矩阵；w_n,1为第n个频率通道的第1个快拍数据对应的系数向量；w_n,l为第n个频率通道的第l个快拍数据对应的系数向量；w_n,L为第n个频率通道的第L个快拍数据对应的系数向量；

为整个探测角度离散化的网格点数；

步骤3：将第n个频率通道的快拍数据对应的稀疏模型通过多频联合稀疏贝叶斯学习算法估计的目标到达角；

所述步骤3实现包括以下子步骤：

步骤3.1：构建系数矩阵

服从的概率分布模型为：

式中，

为所有频率通道的系数矩阵；W₁为第1个频率通道对应的系数矩阵；W_n为第n个频率通道对应的系数矩阵；W_N为第N个频率通道对应的系数矩阵；Γ＝diag(γ)为协方差矩阵；diag(γ)表示以γ作为对角元素的对角矩阵；

γ_i≥0为方差，又称为超参数先验，控制着模型的稀疏性；

表示复高斯分布符号。MFJSBL算法假定所有系数矩阵共享超参数等同于为系数矩阵放置相同的稀疏约束，以增加DOA的拟合概率。

步骤3.2：估计超参数先验γ；

在稀疏贝叶斯学习框架下，依据第Ⅱ类最大似然法，通过使证据函数

相对于γ等于零，并引入因子

得到第i个超参数的更新迭代方程为：

式中，

和

表示所有频率通道的快拍数据；Y₁表示第1个频率通道的快拍数据；Y_n表示第n个频率通道的快拍数据；Y_N表示第N个频率通道的快拍数据；β表示所有频率通道的噪声方差；

表示第1个频率通道的噪声方差；

表示第n个频率通道的噪声方差；

表示第N个频率通道的噪声方差；R_n为第n个频率通道的采样协方差矩阵；

为第n个频率通道的协方差矩阵；

为原先的第i个超参数；

为更新的第i个超参数；a_n,i为第n个频率通道的第i个网格对应的方向系数；上标H是共轭转置运算符，tr(·)表示矩阵的迹；

步骤3.3：估计噪声方差β；

同样依据第Ⅱ类最大似然法，通过使证据函数

的导数相对于β的导数等于零得到第n个频率通道的噪声方差

的更新方程为

式中，

为由Φ_n中与γ中η个最大的峰对应的原子向量组成的索引集

构成。如果目标数K是已知的，峰的数目η＝K；否则，估计的目标数依据Akaike信息准则得到。

步骤3.4：计算算法收敛率；

算法收敛率∈定义为超参数γ的相对改变率，写为

当算法收敛率达到预设的值时停止计算。

步骤3.5：依据得到的超参数γ得到估计的目标到达角；

算法实现时，需首先等间隔离散整个探测区域为

的角度网格。超参数γ值离散的角度网格相对应，因此，依据估计的γ值可得到目标的到达角；

超参数γ中的非零项的网格索引对应的网格角度为目标到达角，而超参数γ中的零元素项表示该网格角没有目标；具体实施时，由于噪声的存在，除了目标所对应的项有非零元素外，其他位置项也可能存在一些小的非零元素。在目标数没有已知的条件下，可依据Akaike信息准则得到估计的目标数。

根据上述实施的具体步骤得到本发明的仿真结果。图3为监测天线是间隔1.5米的ULA，快拍数为L＝15，仿真场景中添加表1中的目标1，当频率通道数n＝1,3,5时，DOA估计的RMSE与目标SNR的关系。可以看出，在不同的选择频带数下，DOA估计的RMSE都随着目标SNR的增加而降低。图4为监测天线是间隔1.5米的ULA，快拍数为L＝15，仿真场景中添加表1中的目标1，目标2或目标1，目标3时，DOA估计的RMSE与目标SNR的关系。可以看出，两个目标估计的RMSE与目标SNR的关系，具有与单目标相同的变化规律。但两个近间隔目标的估计精度要低于两个远间隔目标的。同时，仿真测试了快拍数对估计精度的影响。图5为监测天线是间隔1.5米的ULA，单快拍下，仿真场景中添加表1中的目标1，当频率通道数n＝1,3,5时，DOA估计的RMSE与目标SNR的关系。可以看出，要得到与快拍数L＝15近似相同的RMSE，单快拍下目标SNR需提高15dB左右。为测试阵列形式对目标DOA的影响，图6给出同样阵元数的半径为1.5米的UCA，单快拍下，仿真场景中添加表1中的目标1，当频率通道数n＝1,3,5时，DOA估计的RMSE与目标SNR的关系。可以看出，要达到与ULA(即图5)近似相同的估计精度，目标1的SNR需再提高近似10dB。需指出的是，这里的目标SNR都指的都是相干积累后的目标SNR。总的来说，在不同的目标数、快拍数和阵列形式下，估计精度都随着选择频带数的增多而提高，验证了MFJSBL算法的正确性。

表1：仿真参数

备注：CNR是杂噪比的缩写。

图7为运行1000次蒙特卡洛仿真后，DOA估计的RMSE与目标SNR的关系，其中，频率通道数n＝3和快拍数L＝15。可见，MFJSBL算法的RMSE是五种方法中最低的，更接近于克拉美罗界(CRB)。

图8为在积累的目标SNR为15dB，频率通道数为n＝3和快拍数为L＝15条件下，运行蒙特卡洛仿真1000次得到的DOA估计的RMSE与两目标角度间隔的关系。可以看出，MFJSBL算法在各个角度间隔下都有最好的估计精度。平均的基于期望最大的稀疏贝叶斯学习(AEMSBL)算法对两个目标的估计精度仅次于MFJSBL算法，但当两个目标的夹角间隔为5°时，其性能不稳定。聚焦的基于期望最大的稀疏贝叶斯学习(FEMSBL)和聚焦的基于后验概率最大的稀疏贝叶斯学习(FMAPSBL)算法的性能不如AEMSBL算法。

图9为当数据相干积累1秒，提出算法使用单快拍数据估计的目标DOAs的归一化功率谱。仿真时，假定表1中三个目标在同一距离多普勒元，即目标位于双基距离100km和多普勒频率50Hz处。为了紧凑，仿真调整目标3的DOA为θ＝-20°。可以看出，常规的MUSIC算法在单快拍下是失败的，而MFJSBL算法可正常运行，且估计的目标DOAs与红圈表示三个目标的真实值一致。结果验证了MFJSBL算法在单快拍下估计多目标的可行性。

根据上述实施的具体步骤得到本发明的实测数据分析结果。具体实验于2018年8月在中国德阳开展了飞机目标真实数据的采集工作。实验中，参考天线采用高方向性的八木天线，监测天线采用半径为1.5米的半波偶极子天线组成的7元UCA。采用基带信号采样率为200KHz的数字滤波器，从接收的宽带数据中提取成都猛追湾电视塔的97MHz和102.6MHz的通道数据。基带信号采样率设置为200KHz，数据积累时间设为2秒。民航飞机DOA的真实值由测量活动中获得的自动相关监测广播(ADS-B)记录确定。

表2：实验数据估计的目标DOA的RMSE单位：(°)

表2总结了目标DOA估计的RMSE统计结果。结果显示，在DOA估计误差方面，MFJSBL算法性能优于与使用单个频率通道数据的MAPSBL算法和联合两个频率通道数据的AMAPSBL算法。这说明MFJSBL算法可在一定程度上有效提高目标的估计精度。对此有两种可能的解释：一是联合两个频率通道的数据在一定程度上可提高目标的等效SNR。二是MFJSBL采用约束超参数策略拟合DOA，使算法性能即使在小干扰下也表现稳定。AMAPSBL算法性能可能会受到频率通道数据中某些特定时刻差的估计性能影响比较大。另外，实际数据的DOA估计结果相对较差是因为实际数据的目标SNR与达到图6所示的估计精度所需的目标SNR相比仍是较低的。总体来讲，实际数据结果与仿真分析相吻合，验证了MFJSBL算法的正确性。

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

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

Claims (2)

1.一种基于多频联合稀疏贝叶斯学习的到达角估计方法，其特征在于，具体包含以下步骤：

步骤1：获取多个频率通道的快拍数据；

步骤2：构建步骤1所述的第n个频率通道快拍数据对应的稀疏模型，构建利用所有频率通道的快拍数据联合估计到达角的稀疏模型；

步骤3：将第n个频率通道的快拍数据对应的稀疏模型通过多频联合稀疏贝叶斯学习算法估计的目标到达角；

步骤2所述构建第n个频率通道快拍数据对应的基于网格的稀疏模型为：

Y_n＝Φ_nW_n+Z_n n＝1,...,N

式中，

W_n＝[w_n,1,...,w_n,l,...,w_n,L]；Φ_n为第n个频率通道的过完备字典矩阵；

为第1个网格角

对应的第n个频率通道的阵列流型向量；

为第k个网格角

对应的第n个频率通道的阵列流型向量；

为第

个网格角

对应的第n个频率通道的阵列流型向量；W_n为第n个频率通道的过完备字典矩阵对应的所有快拍数据的系数矩阵；w_n,1为第n个频率通道的第1个快拍数据对应的系数向量；w_n,l为第n个频率通道的第l个快拍数据对应的系数向量；w_n,L为第n个频率通道的第L个快拍数据对应的系数向量；

为整个探测角度离散化的网格点数；

所述步骤3实现包括以下子步骤：

步骤3.1：构建系数矩阵

服从的概率分布模型为：

式中，

为所有频率通道的系数矩阵；W₁为第1个频率通道对应的系数矩阵；W_n为第n个频率通道对应的系数矩阵；W_N为第N个频率通道对应的系数矩阵；Γ＝diag(γ)为协方差矩阵；diag(γ)表示以γ作为对角元素的对角矩阵；

为方差，又称为超参数先验，控制着模型的稀疏性；

表示复高斯分布符号；MFJSBL算法假定所有系数矩阵共享超参数等同于为系数矩阵放置相同的稀疏约束，以增加DOA的拟合概率；

步骤3.2：估计超参数先验γ；

在稀疏贝叶斯学习框架下，依据第Ⅱ类最大似然法，通过使证据函数

相对于γ等于零，并引入因子

得到第i个超参数的更新迭代方程为：

式中，

和

表示所有频率通道的快拍数据；Y₁表示第1个频率通道的快拍数据；Y_n表示第n个频率通道的快拍数据；Y_N表示第N个频率通道的快拍数据；β表示所有频率通道的噪声方差；

表示第1个频率通道的噪声方差；

表示第n个频率通道的噪声方差；

表示第N个频率通道的噪声方差；R_n为第n个频率通道的采样协方差矩阵；

为第n个频率通道的协方差矩阵；

为原先的第i个超参数；

为更新的第i个超参数；a_n,i为第n个频率通道的第i个网格对应的方向系数；上标H是共轭转置运算符，tr(·)表示矩阵的迹；

步骤3.3：估计噪声方差β；

同样依据第Ⅱ类最大似然法，通过使证据函数

的导数相对于β的导数等于零得到第n个频率通道的噪声方差

的更新方程为

式中，

为由Φ_n中与γ中η个最大的峰对应的原子向量组成的索引集

构成；如果目标数K是已知的，峰的数目η＝K；否则，估计的目标数依据Akaike信息准则得到；

步骤3.4：计算算法收敛率；

算法收敛率∈定义为超参数γ的相对改变率，写为

当算法收敛率达到预设的值时停止计算；

步骤3.5：依据得到的超参数γ得到估计的目标到达角；

算法实现时，需首先等间隔离散整个探测区域为

的角度网格；超参数γ值离散的角度网格相对应，因此，依据估计的γ值可得到目标的到达角。

2.根据权利要求1所述的基于多频联合稀疏贝叶斯学习的到达角估计方法，其特征在于，所述步骤1具体实现包括以下子步骤：

步骤1.1：构建第t个时刻采集的第m个监测通道接收到的第n个频率通道的信号，具体如下：

式中，

表示第t个时刻采集的第n个频率通道的参考信号；

为直达波信号相对于参考信号的归一化幅度；

为直达波相对于参考信号的时延；

为直达波的到达角；

表示第p条多径信号相对于参考信号的复幅度；

表示第p条多径信号相对于参考信号的时延；

表示第p条多径信号的DOA；η_k表示第k个目标的幅度；τ_k表示第k个目标的时延；

分别表示第k个目标的多普勒频率；θ_k表示第k个目标的DOA；

表示第n个频率通道接收的第m个监测信号中的高斯白噪声；

表示监测信号中第n个频率通道的目标、直达波和多径信号的阵列导向矢量；M表示监测通道数，N表示频率通道数，K表示目标数，N_c表示多径数，T_int＝N_intT_s表示每场数据积累时间，N_int表示积累数据点数，T_s＝1/F_s表示离散时间间隔，F_s表示采样率；

步骤1.2：将第t个时刻采集的第m个监测通道接收到的第n个频率通道的信号利用杂波抑制算法滤除监测信号中明显强于目标的直达波和多径杂波，得到抑制后第t个时刻第m个监测通道接收的第n个频率通道的信号为：

步骤1.3：通过第m个监测通道接收的第n个频率通道的信号与相应的参考信号互相关生成的距离多普勒谱为：

τ＝[0,T_s,...,(R_max/ΔR)*T_s]；f_d＝[(-v_max/Δv)*Δf,...,0,...,(v_max/Δv)*Δf]

式中，(·)^*为共轭运算符；R_max表示外源雷达可探测的最大距离；v_max表示外源雷达可探测的最大速度；ΔR＝c/B表示可分辨的距离间隔，c表示光速，B表示照射源信号带宽；Δv＝λ*Δf表示可分辨的速度间隔，λ表示照射源信号波长，Δf＝1/T_int表示可分辨的多普勒间隔；

步骤1.4：利用恒虚警率检测方法获取第m个监测通道接收的第n个频率通道的信号与相应的参考信号互相关生成的距离多普勒谱中目标的复幅度值作为到达角估计的快拍数据；假定K个目标位于同一距离多普勒，则第m个监测通道接收的第n个频率通道的信号与相应的参考信号互相关生成的距离多普勒谱中第k个时延τ_k和第k个多普勒

处的目标幅度值为：

式中，

表示第n个频率通道接收的第k个目标在时延和多普勒

处的自模糊函数值；

表示第n个频率通道接收的第m个监测信号相干积累后的噪声；

步骤1.5：整理估计目标DOA所需的快拍数据；这里需指出的，步骤1.4直接获取的第n个频率通道的快拍数据为

是单个的快拍数据；

其中，

为增加估计精度，采用频域滑窗的方式计算步骤1.4的距离多普勒谱L次，得到L个快拍数据；

整理第n个频率通道的L个快拍数据为矩阵形式为：

Y_n＝A_nC_n+Z_n n＝1,...,N

式中，

Y_n＝[y^n,1,...,y^n,l,...,y^n,L]，

A_n＝[a_n(θ₁),...,a_n(θ_k),...,a_n(θ_K)]，

Y_n表示第n个频率通道对应的所有监测通道的快拍数据；y^n,1表示第n个频率通道的第1快拍在所有监测通道上的数据；y^n,l表示表示第n个频率通道的第l快拍在所有监测通道上的数据；y^n,L表示表示第n个频率通道的第L快拍在所有监测通道上的数据；

表示第n个频率通道的第1个监测通道上获取的第l个快拍数据；

表示第n个频率通道的第m个监测通道上获取的第l个快拍数据；

表示第n个频率通道的第M个监测通道上获取的第l个快拍数据；上标T为转置运算符，L为快拍数；A_n表示第n个频率通道接收的所有目标对应的阵列流型矩阵；a_n(θ₁)表示第n个频率通道接收的第1个目标对应的阵列流型矢量；a_n(θ_k)表示第n个频率通道接收的第k个目标对应的阵列流型矢量；a_n(θ_K)表示第n个频率通道接收的第K个目标对应的阵列流型矢量；

表示第n个频率通道接收的第k个目标在第1个监测通道的阵列流型；

表示第n个频率通道接收的第k个目标在第m个监测通道的阵列流型；

表示第n个频率通道接收的第k个目标在第M个监测通道的阵列流型；C_n为第n个频率通道接收目标的复幅度；Z_n为第n个频率通道的高斯白噪声。

Images