CN110806557A - 一种冲击噪声背景下单基地mimo雷达doa估计方法 - Google Patents

Abstract

本发明公开一种冲击噪声背景下单基地MIMO雷达DOA估计方法，第一步，基于SαS稳态分布的冲击噪声建模以及基于柯西核函数的冲击噪声抑制；第二步，基于迭代凸优化算法进行信号子空间构建；第三步，基于旋转不变方法来实现DOA估计。本发明针对Cauchy核函数的求解问题，利用复值牛顿梯度下降算法实现了对阵列接收数据信号子空间的迭代求解；经过多次迭代操作，得到性能良好的信号子空间。利用Cauchy核函数实现了对冲击噪声的抑制，解决了单基地MIMO雷达中冲击噪声背景下DOA估计问题。方位参数信息估计结果较为精确；在避免谱峰搜索的同时实现良好的来波方位估计。

Description

一种冲击噪声背景下单基地MIMO雷达DOA估计方法

【技术领域】

本发明一种冲击噪声背景下单基地MIMO雷达DOA估计方法，属于雷达探测技术领域，具体涉及单基地多输入多输出(MIMO)雷达的无源定位，特别是在冲击噪声背景下的单基地MIMO雷达的来波方位(DOA)估计技术。

【背景技术】

MIMO雷达是一种新体制雷达，其通过在发射端发射相互正交的信号，并在接收端进行正交匹配滤波将正交信号分离，能够在空间形成多个观测通道，且能够抑制环境噪声、杂波等非理性信号，在自身生存能力和目标探测能力方面都有很大的提高，具有广阔的应用前景。

目前，对MIMO雷达的DOA估计研究主要基于高斯噪声假设而展开，但在很多应用场景下常规高斯噪声假设实际上难以成立。如海杂波、地杂波和逆向散射回波对DOA估计的影响常表现出冲击特性，且服从冲击噪声分布的随机变量在理论上不具有二阶以上的统计量。因此，常规的基于二阶相关函数和高阶累积量的信号处理方法在冲击噪声背景下性能严重退化，甚至失效。在上述应用场景中如何有效地抑制冲击噪声对MIMO雷达DOA估计的影响，面临着两个主要的技术难题：一是针对冲击噪声的冲击特性，实现对冲击噪声中异常值抑制；二是异常值抑制后，需要高效算法实现冲击噪声背景下的MIMO雷达DOA精确估计。

【发明内容】

本发明一种冲击噪声背景下单基地MIMO雷达DOA估计方法，要解决的技术问题是：针对单基地MIMO雷达，利用柯西核函数实现对单基地MIMO雷达冲击噪声抑制，结合迭代凸优化算法和信号子空间旋转不变方法完成单基地MIMO雷达DOA精确估计。

本发明采取的技术方案如下：

第一步，基于SαS稳态分布的冲击噪声建模以及基于柯西(Cauchy)核函数的冲击噪声抑制

首先针对冲击噪声的建模不唯一性，这里采用SαS稳态分布来构建冲击噪声模型，并且为柯西核函数估计器对冲击噪声的抑制提供相应的数据模型。假设冲击噪声服从SαS稳态分布，SαS稳态分布的特征函数

可以表示为

其中

j是单位虚数，exp是指数函数，α表示特征指数，其控制SαS稳态分布的厚度，γ表示散度参数，其类似于高斯分布的方差，β表示对称参数，β＝0意味着对称SαS稳态分布，a表示位置参数，意味着分布的均值或者中值。越小的α，冲击噪声的强度越强。

标准SαS稳态分布的概率密度函数可以表示为

其中

不同特征指数的核函数可以定义为

利用柯西核函数来实现对冲击噪声的抑制，基于柯西核函数的DOA估计器可表示为

其中M表示单基地MIMO雷达发射阵元个数，N表示接收阵元个数，argmin表示最小化，e_i(1≤i≤M+N-1)表示系统误差，ξ表示待估计参数。

第二步，基于迭代凸优化算法进行信号子空间构建

为了实现对构建的柯西核函数DOA估计器进行精确的求解，这里采用复值牛顿梯度下降迭代凸优化算法来实现对原始非凸优化问题的求解。首先把原始的信号模型分解成相应的信号子空间和噪声子空间，然后采用交替迭代牛顿梯度下降算法对信号子空间和噪声子空间进行精确重构。具体步骤如下：

首先，对单基地MIMO的阵列接收信号Y进行低秩分解

Y＝UV (公式六)

其中U是(M+N-1)×Q的列满秩矩阵，V是Q×L的行满秩矩阵。在冲击噪声背景下的剩余残差矩阵可以表示为Y-UV＝[r₁,r₂,...,r_l,...,r_L] (公式七)

其中r_l＝y_l-Uv_l表示剩余残差，y_l和v_l分别表示Y和V的第l列。

然后，为实现U和V的求解，定义以下目标函数J_ρ

其中ρ(·)表示柯西核M估计器，min表示最小化操作，[q₁,...,q_M+N-1]＝(Y-UV)^T＝(Y^T-V^TU^T)表示(Y-UV)^T的相应的列矢量。进而采用如下迭代凸优化算法

其中V^(k+1)和U^(k+1)分别表示第k+1的迭代结果。

从以上等式可以发现，对于U和V的迭代求解具有相同的处理过程，为简化分析，这里只讨论对于V的求解。对于V的迭代凸优化求解过程可以划分为相应的L个子问题

其中，

表示

的第l列，令r^(k+1)＝y_l-U^(k)v_l，则相应的柯西核代价函数可表示为

为实现对的求解，采用如下的复值牛顿算法

其中，

W是对角阵

因此，f关于复矢量

的Hessian矩阵可以表示为

其中，

表示关于r^(k+1)的(M+N-1)×(M+N-1)Hessian矩阵

可以通过以下的迭代得到

其中

并且步长μ_k可通过黄金分割搜索法得到

当U和V的迭代更新满足如下的条件时，迭代终止

第三步，基于旋转不变方法来实现DOA估计

通过迭代凸优化算法的多次迭代，得到的维度为(M+N-1)×Q的U矩阵对应于阵列接收数据相应的信号子空间。在得到信号子空间以后，旋转不变技术可以实现高效精确的方位参数求解。选择包含M+N-2个阵元的一对子阵Sub_X＝{0,1,2,…,M+N-3}和Sub_Y＝{1,2,3,…,M+N-2}，则两个子阵相应的信号子空间可以表示为

其中U_X和具有相同的信号子空间，同理U_Y和

也具有相应的信号子空间。因此，一个Q×Q维矩阵F满足如下的关系

由于U_X和U_Y具有相同的信号子空间，增广矩阵U_XY＝[U_X U_Y]和U_X，U_Y具有相同的秩。因此，存在一个2Q×Q维正交矩阵P＝[P_X P_Y]^T满足：

通过以上分析可以发现U_X和U_Y满足如下的关系

U_Y＝U_XΨ (公式二十四)

其中，

选择矩阵Φ和算子Ψ满足如下关系

Φ＝FΨF^-1 (公式二十五)

从公式可以看出，Φ和算子Ψ具有相同的奇异值，所以最终的方位参数可以直接从算子Ψ中得到

其中ψ_q表示Ψ的第q个奇异值，imag表示复数的虚部。

本发明的有益效果主要包括：

第一，针对Cauchy核函数的求解问题，利用复值牛顿梯度下降算法实现了对阵列接收数据信号子空间的迭代求解。经过多次的迭代操作，能够得到性能良好的信号子空间。

第二，利用Cauchy核函数实现了对冲击噪声的抑制，解决了单基地MIMO雷达中冲击噪声背景下的DOA估计问题。

第三，方位参数信息估计结果较为精确。在得到信号子空间之后，旋转不变技术被用来实现对方位参数的求解。所提出的旋转不变技术在避免谱峰搜索的同时实现良好的来波方位估计。

【附图说明】

图1a～图1d是不同特征指数的冲击噪声示意图。

图2是SαS分布核函数示意图。

图3是Cauchy核M估计器随误差的变化示意图。

图4是SC-ESPRIT空间谱示意图。

图5是RMSE随快拍数变化的示意图。

图6是RMSE随GSNR变化的示意图。

图7是RMSE随角度间隔变化的示意图。

图8是本发明方法流程框图。

【具体实施方式】

下面结合附图对本发明作进一步的说明。本发明所使用的软件环境是Matlab R2017a版本,所使用的硬件计算环境是Intel Corel i7处理器，主频2.4GHz，内存8GHz。仿真中单基地MIMO雷达的发射天线M＝5，接收天线N＝6，发射和接收天线均是阵元间距为半个波长的均匀线阵。稳态分布α的变化范围为0＜α≤2，由于0＜α≤1在实际情况中很难碰到，所以本发明只考虑1＜α≤2情况下的冲击噪声。由于冲击噪声不存在有限的二阶矩，所以仿真中的信噪比均采用广义信噪比(GSNR)的形式且GSNR的变化范围为0dB～10dB。快拍数(snapshots)的变化范围是20～200。在军方误差(RMSE)随GSNR的变化的算法性能对比时，快拍数的选取固定为100。在RMSE随snapshots的变化的算法性能对比时，GSNR的选取固定为5dB。在角度分辨率性能的试验验证中，角度间隔Δ的变化范围是1°～15°。且此时GSNR和snapshots分别固定为5dB和100。其它具体的仿真参数在各个子实验中分别进行了详细的说明。如图8所示，本发明方法步骤如下：

第一步，基于SαS稳态分布的不同特征指数α下的冲击噪声示意图以及柯西核DOA估计器的误差分析。如图1a～图1d所示，冲击噪声相对于高斯噪声具有较大的的幅度。并且越小的α，冲击噪声的强度越强。图2可以看出，对于不同特征指数的稳态分布，其核函数呈现出类似的特征，即，在随机变量的幅度较小时，核函数单调递增，随着随机变量幅度的增加，核函数逐渐趋近于1/x。由于α＝1对应的柯西分布核函数具有简洁的表达式，所以这里采用柯西核函数来实现对噪声的抑制。并且从图3可以看出随着误差的增大，柯西核函数估计器的变化相对比较平缓，这也说明了柯西核函数有助于实现对噪声的抑制。

第二步，利用迭代凸优化算法对柯西核函数估计器进行多次迭代求解，最终得到相应的信号子空间。

第三步，然后再基于旋转不变理论实现对单基地MIMO雷达角度参数的精确估计。如图4所示，仿真中入射的窄带非相关信源数K＝4，且入射方向来为[-30°,-10°,10°,30°]，特征指数α取1.5来说明相对强的冲击噪声的影响。GSNR设置为5dB，快拍数取100。从图4中可以看出，在较强的冲击噪声的背景下，所提出的算法仍能够实现良好的方位估计。

进一步的，以下为冲击噪声背景下单基地MIMO雷达不同DOA估计算法RMSE性能对比。如图5所示，首先考虑RMSE性能随快拍数的变化，通过对比发现所提出的SC-ESPRIT算法在不同的快拍数下具有较低的误差。然后考虑RMSE性能随GSNR的变化，如图6所示，随着广义信噪比的增加，所提出的SC-ESPRIT算法性能逐渐趋近CRB。最后，考虑RMSE性能随角度间隔的变化。如图7所示，在角度间隔较小时(1≤Δθ≤4°)，所提出的SC-ESPRIT算法仍然能够实现良好的方位估计。当角度间隔较大时(7≤Δθ≤15°)，所提出的SC-ESPRIT算法性能仍然优越于其他算法。通过图5-图7的RMSE对比可以发现，所提出的SC-ESPRIT算法在快拍较少、低信噪比和信源间隔较小的情况下都能实现良好的估计性能。

Claims (4)

1.一种冲击噪声背景下单基地MIMO雷达DOA估计方法，其特征在于：该方法步骤如下：

第一步，基于SαS稳态分布的冲击噪声建模以及基于柯西核函数的冲击噪声抑制

首先针对冲击噪声的建模不唯一性，采用SαS稳态分布来构建冲击噪声模型，并且为柯西核函数估计器对冲击噪声的抑制提供相应的数据模型；

第二步，基于迭代凸优化算法进行信号子空间构建

为了实现对构建的柯西核函数DOA估计器进行精确的求解，采用复值牛顿梯度下降迭代凸优化算法来实现对原始非凸优化问题的求解：首先把原始的信号模型分解成相应的信号子空间和噪声子空间，然后采用交替迭代牛顿梯度下降算法对信号子空间和噪声子空间进行精确重构；

第三步，基于旋转不变方法来实现DOA估计

通过迭代凸优化算法的多次迭代，得到的维度为(M+N-1)×Q的U矩阵对应于阵列接收数据相应的信号子空间；在得到信号子空间以后，通过旋转不变技术实现高效精确的方位参数求解。

2.根据权利要求1所述的一种冲击噪声背景下单基地MIMO雷达DOA估计方法，其特征在于：所述第一步的具体过程如下：

假设冲击噪声服从SαS稳态分布，SαS稳态分布的特征函数

可以表示为

其中

j是单位虚数，exp是指数函数，α表示特征指数，其控制SαS稳态分布的厚度，γ表示散度参数，其类似于高斯分布的方差，β表示对称参数，β＝0意味着对称SαS稳态分布，a表示位置参数，意味着分布的均值或者中值；越小的α，冲击噪声的强度越强；

标准SαS稳态分布的概率密度函数可以表示为

其中

不同特征指数的核函数可以定义为

利用柯西核函数来实现对冲击噪声的抑制，基于柯西核函数的DOA估计器可表示为

其中M表示单基地MIMO雷达发射阵元个数，N表示接收阵元个数，argmin表示最小化，e_i(1≤i≤M+N-1)表示系统误差，ξ表示待估计参数。

3.根据权利要求1所述的一种冲击噪声背景下单基地MIMO雷达DOA估计方法，其特征在于：所述第二步的具体过程如下：

首先，对单基地MIMO的阵列接收信号Y进行低秩分解

Y＝UV (公式六)

其中U是(M+N-1)×Q的列满秩矩阵，V是Q×L的行满秩矩阵；在冲击噪声背景下的剩余残差矩阵可以表示为

Y-UV＝[r₁,r₂,...,r_l,...,r_L] (公式七)

其中r_l＝y_l-Uv_l表示剩余残差，y_l和v_l分别表示Y和V的第l列；

然后，为实现U和V的求解，定义以下目标函数J_ρ

其中ρ(·)表示柯西核M估计2器，min表示最小化操作，[q₁,...,q_M+N-1]＝(Y-UV)^T＝(Y^T-V^TU^T)表示(Y-UV)^T的相应的列矢量；进而采用如下迭代凸优化算法

其中V^(k+1)和U^(k+1)分别表示第k+1的迭代结果；

从以上等式可以发现，对于U和V的迭代求解具有相同的处理过程，为简化分析，这里只讨论对于V的求解；对于V的迭代凸优化求解过程可以划分为相应的L个子问题

其中，

表示

的第l列，令r^(k+1)＝y_l-U^(k)v_l，则相应的柯西核代价函数可表示为

为实现对

的求解，采用如下的复值牛顿算法

其中，

W是对角阵

因此，f关于复矢量

的Hessian矩阵可以表示为

其中，

表示关于r^(k+1)的(M+N-1)×(M+N-1)Hessian矩阵

可以通过以下的迭代得到

其中

并且步长μ_k可通过黄金分割搜索法得到

当U和V的迭代更新满足如下的条件时，迭代终止

4.根据权利要求1所述的一种冲击噪声背景下单基地MIMO雷达DOA估计方法，其特征在于：所述第三步的具体过程如下：

选择包含M+N-2个阵元的一对子阵Sub_X＝{0,1,2,,M+N-3}和Sub_Y＝{1,2,3,,M+N-2}，则两个子阵相应的信号子空间可以表示为

其中U_X和A_X＝A(1:M+N-2,:)具有相同的信号子空间，同理U_Y和A_Y＝A(2:M+N-1,:)也具有相应的信号子空间；因此，一个Q×Q维矩阵F满足如下的关系

U_X＝A_XF (公式二十一)

U_Y＝A_YF＝A_XΦF (公式二十二)

由于U_X和U_Y具有相同的信号子空间，增广矩阵U_XY＝[U_X U_Y]和U_X，U_Y具有相同的秩；因此，存在一个2Q×Q维正交矩阵P＝[P_X P_Y]^T满足：

通过以上分析可以发现U_X和U_Y满足如下的关系

U_Y＝U_XΨ (公式二十四)

其中，选择矩阵Φ和算子Ψ满足如下关系

Φ＝FΨF^-1 (公式二十五)

从公式可以看出，Φ和算子Ψ具有相同的奇异值，所以最终的方位参数可以直接从算子Ψ中得到

其中ψ_q表示Ψ的第q个奇异值，imag表示复数的虚部。

Images