一种稳健的大阵列MIMO雷达目标收发角联合估计方法与流程

本发明属于多输入多输出(mimo)雷达目标参数估计领域；涉及收发阵元数与快拍数相接近的大阵列mimo雷达系统；涉及非高斯噪声背景下双基地mimo雷达的目标收发角联合估计；涉及收发阵元数大于快拍数情况下的mimo雷达目标收发角度估计。

背景技术：

雷达(radar)通过发射无线电波来检测目标并测定目标的空间位置，从而对目标进行定位。作为一种新体制雷达，多输入多输出(mimo)雷达由多个发射天线同时发射多个信号探测目标，并使用多个接收天线接收目标回波信号，通过采用多种可区分的发射信号(本发明采用正交信号)，为雷达引入参数上有差异的目标观测通道。由于各收发通道之间采用相干处理，mimo雷达可充分利用收发阵元获得虚拟孔径，提高收发角估计的性能。

mimo雷达的目标收发角联合估计问题，是指在噪声背景下对阵列接收的回波信号进行有效处理，进而对目标的发射角(dod)和接收角(doa)进行估计。根据发射接收天线的阵列布局不同，mimo雷达可分为分布式mimo雷达和集中式mimo雷达两类。分布式mimo雷达的发射天线之间距离较远以获得空间分集增益；集中式mimo雷达的发射阵列和接收阵列的各个阵元间距通常在半波长以内，通过相干处理提高估计增益，比如单/双基地mimo雷达。在双基地mimo雷达中，利用收发角联合估计可实现目标定位。现有的收发角联合估计方法主要有2d-capon方法、2d-music方法、分维esprit方法等。这些方法虽然能有效估计出收发角度，但仍存在很多问题亟需解决，主要包括：(1)现有方法通常将噪声模型简化为理想的高斯噪声，其分布函数为正态分布；而在实际环境中，噪声不是简单的高斯噪声，而是表现为一些非高斯特性，比如服从重尾分布，若仍采用高斯噪声背景下的目标参数估计方法，将导致估计性能严重下降；(2)现有方法的前提条件都是假定收发阵元数远远小于快拍数，而当收发阵元数大于快拍数时，在求解采样协方差矩阵定点迭代的过程中会出现因秩缺失而无解的问题，此时无法对收发角度进行有效估计；(3)现有方法所考虑的mimo雷达系统均假设快拍数足够多，而实际中，当快拍数与收发阵元数相接近时，现有估计方法不再适用。在实际的雷达例如超视距mimo雷达系统中，往往会出现收发阵元数与快拍数相接近甚至收发阵元数大于快拍数的情况，即大阵列系统。迄今为止，非高斯噪声背景下大阵列mimo雷达系统的目标收发角联合问题一直是一个难题。

近年来，随机矩阵理论在金融数学、核物理、数字通信等领域应用日益广泛，已成为学术界的一个新的研究热点。然而，随机矩阵理论在大阵列mimo雷达目标参数估计领域的应用研究甚少，尤其是在非高斯噪声背景且收发阵元数大于快拍数情况下，针对大阵列双基地mimo雷达的目标收发角联合估计方法，国内外尚未见相关成果的研究报道。

随机矩阵渐进谱理论是符合理想高斯噪声背景下大阵列mimo雷达目标收发角联合估计的有效方法。然而，实际环境中，噪声往往不再服从正态分布，而是呈现重尾分布。针对非高斯噪声情况，由于采样协方差矩阵将不再是统计协方差矩阵的一致性估计，导致高斯噪声背景下提出的收发角联合估计方法不再适用。为解决这一问题，人们利用极大似然估计对采样协方差矩阵进行定点迭代，作为协方差矩阵的一致估计。然而，当阵元数大于快拍数时，定点迭代的过程中会出现因秩缺失而无解的情况，这一问题是目前亟待解决的难题。

近年来出现的基于线性收缩(linearshrinkage)的协方差矩阵估计方法，已成为一种新的信号处理手段。收缩估计可以在矩阵维数大于快拍数的情况下，对大维采样协方差矩阵的特征值分布区间进行收缩，解决秩缺失的问题。目前，已经有文献研究了大维系统中协方差矩阵的线性收缩估计问题以及一些改进的方法，广义收缩估计器已经应用于空时自适应处理(stap)检测问题中。然而，在大阵列mimo雷达收发角联合估计领域，国内外尚未见任何有关线性收缩应用的报道。

技术实现要素：

本发明提供一种稳健的大阵列mimo雷达目标收发角联合估计方法，以解决估计性能严重下降、无法对收发角度进行有效估计、非高斯噪声背景下大阵列mimo雷达系统的目标收发角联合估计困难的问题。

本发明采取的技术方案是：建立双基地mimo雷达的系统模型，m个发射阵元同时发射窄带信号照射k个远场目标，在接收端采用n个接收阵元，其中发射阵列和接收阵列均采用均匀线阵，阵元间距dt＝dr＝λ/2，发射角度为θ，接收角度为第k个目标的发射角度和接收角度分别为θk和l为快拍数，其特征在于包括下列步骤：

步骤一、mimo雷达接收端n个接收阵元对多次快拍得到的回波数据信号进行收集和处理，第l个快拍期间回波数据可以写成一个n×lt的矩阵x^(l)：

其中，表示m行1列的发射方向矢量，t表示发射端，表示n行1列的接收方向矢量，r表示接收端；s表示发射波形矩阵，这里sm(m＝1,…,m)表示第m个发射阵元所发射的波形，其编码长度为lt，且s为归

一化正交矩阵，ss^h＝im；表示在第l个快拍期间第k个目标的散射系数；

表示一个n行lt列的非高斯噪声矩阵，服从student-t分布；

步骤二、首先对n个接收阵元所接收的信号x^(l)进行匹配滤波处理，得到：

其中，为n行m列的匹配滤波后的信号；

步骤三、对匹配滤波后的结果进行矢量化处理，即y^(l)＝vec(x^(l)s^h)，在第l个快拍期间，经过预处理后的观测数据可表示为：

其中，为目标散射系数矩阵；是在第l个快拍期间处理后的噪声矢量，维数为mn行1列；是联合导向矢量矩阵，维数为mn行k列at(θk)和分别为发射和接收方向矢量，是khatri-rao矩阵积，是kronecker矩阵积；

步骤四、收集l个快拍期间的数据y^(l)(l＝1,…,l)，可得到一个观测数据矩阵y，它是一个维数为mn×l的随机矩阵，即：

其中，y＝[y⁽¹⁾,…,y^(l)]，b＝[b⁽¹⁾,…,b^(l)]，n＝[n⁽¹⁾,…,n^(l)]，这里表示包含收发角度的联合导向矩阵，b是服从均值为0，方差为1的目标散射系数矩阵，n为mn×l维的服从student-t分布非高斯噪声矩阵，由高斯噪声矩阵的元素乘上服从伽马分布γ(v,1/v)的因子构成，其中，v＝0.02；

步骤五、接收信号的采样协方差矩阵可以表示为求解协方差矩阵的极大似然估计；

步骤六、当收发阵元数大于快拍数时，l<mn，上述极大似然估计方法因为秩缺失而导致无解，此时采用线性收缩方法使得定点迭代的过程能持续进行；

步骤七、在获得稳健协方差矩阵估计的基础上，利用2d-mw方法可以得到收发角联合估计的谱峰搜索函数；

步骤八、当收发阵元数与快拍数在同一数量级时，即mn→∞，l→∞，mn/l＝c，c为一个常数，y可以看作一个大维随机矩阵，传统方法性能急剧下降，不能正确估计出收发角度；因此，应用随机矩阵理论中的渐近谱理论和围线积分技术来产生一个稳健的谱峰搜索函数；

步骤九、对进行谱峰搜索可以获得发射角θ和接收角

本发明所述步骤五的具体实现方法是：对接收信号y^(l)进行归一化yi的概率密度函数为其中u(mn/2)/2π^mn/2为单位圆上的均匀密度分布函数，对p(yi:r)求导并令导数等于零可以得到经过定点迭代可以得到极大似然估计定点迭代的初始矩阵可以是任意正定矩阵，令为单位矩阵。

本发明所述步骤六的收缩方法如下，其中0<ρ<1为收缩系数，将定点迭代完成之后最终收敛的协方差矩阵记为收缩系数ρo由最小均方误差准则得到通过化简可以解得收缩系数为：

由于r未知，用替代r可以得到：

本发明所述步骤七具体实现方法是：对采样协方差矩阵进行特征值分解得到

其中，(λmn>0)，ee^h＝i。定义0<z<∞上的一个连续的函数f(z)，其反函数为可以得到

收发角联合估计估计的谱峰搜索函数可以表示为

改变函数f(z)可以得到不同的谱峰搜索函数，令函数f(z)＝ln(z)，得到midway(mw)方法的代价函数

本发明所述步骤八具体实现方法是：首先根据留数定理，可以表示为：

其中，ζ^-是一条在上的顺时针方向曲线，并且该曲线包含r的所有特征值；

对曲线ζ^-进行参数化处理，令其中，是一个包含所有特征值的矩形，表示该矩形的边界，wmn(z)与特征值有关，其关系式如下：

为了简化公式，定义两个变量bmn(z)和mmn(z)来分别描述特征值和特征向量：

得到参数化表示形式为

通过变量bmn(z)，得到wmn(z)的导数对上式进一步化简得到：

当mn→∞，l→∞，mn/l＝c时，用bmn(z)和mmn(z)来描述特征值和特征向量不再准确，此时需要研究特征值和特征向量的渐近表现，在随机矩阵理论中，采样协方差矩阵特征值的瞬时分布函数的stieltjes变换定义如下：

仅仅体现了采样协方差矩阵的特征值渐近特性，为了同时体现特征值和特征向量的渐近特性，定义如下函数：

此时，我们用和替换bmn(z)和mmn(z)，得到：

其中，

为了求解定义一个关于的函数：

其中，x∈[0,c]；

在x∈(0,c)区间上对求导：

利用留数定理可以得到：

其中，是下面等式的解：

可以由柯西留数定理直接得到，再此基础上对求积分，可以得到如下：

令x＝c，经过化简可以得到ξmn(c)：

其中，是下面等式的解：

至此，基于线性收缩和随机矩阵理论的mimo雷达收发角联合估计谱峰搜索函数可以写为：

其中，如上所述。

本发明在非高斯背景下，且以student-t噪声为例，计算接收信号的采样协方差矩阵，并利用定点迭代技术产生其协方差矩阵的一致估计，即：其中，yi为归一化接收信号向量，j为迭代的次数。当收发阵元数的乘积mn大于快拍数l时，定点迭代的过程中会出现因秩缺失而无法求解采样协方差矩阵的问题，因此本发明采用线性收缩的方法解决这一问题，得到稳健的采样协方差矩阵估计：

其中，ρ为收缩系数，本发明中采用最小均方误差准则对其进行求解。

本发明将基于pisarenko结构的midway(mw)方法应用到mimo雷达收发角联合估计当中，提出一种2d-mw方法，当收发阵元数和快拍数在同一数量级时，采用随机矩阵理论中的围线积分技术、stieltjes变换以及留数定理，推导了大维条件下稳健的代价函数，提出一种基于随机矩阵理论的2d-rmw方法，之后通过谱峰搜索可以得到精确的收发角。在相同仿真条件下，将2d-rmw方法与2d-mw方法进行了对比，验证了方法的有效性，代价函数如下所示：

其中，ξmn(mn＝1,…,mn)为基于随机矩阵理论推导出的谱峰搜索函数的修正值，为联合导向矢量，为特征向量。

本发明在非高斯噪声背景且收发阵元数大于快拍数的情况下，将线性收缩技术与随机矩阵理论结合，通过收缩得到采样协方差矩阵的一致估计，再利用高斯噪声背景下基于随机矩阵理论的2d-rmw方法，提出了一种新的收发角联合估计方法----2d-rlsmw方法，代价函数如下：

其中，为基于线性收缩和随机矩阵理论结合推导出的谱峰搜索函数的修正值，ls指的是线性收缩，为联合导向矢量，为特征向量。

本发明为了分析发明方法的收发角联合估计性能，在大阵列双基地mimo雷达系统中，通过matlab仿真实验证明了高斯环境背景下2d-rmw方法的有效性，检验了student-t分布的非高斯噪声背景下且收发阵元数大于快拍数时本发明所提出的2d-rlsmw方法的优越性。

本发明的优点在于：(1)在大阵列mimo雷达收发角联合估计中，利用极大似然估计，通过定点迭代产生协方差矩阵的一致估计，解决非高斯噪声背景下收发角估计不准确的问题；(2)首次将线性收缩技术应用到大阵列mimo雷达目标收发角联合估计中，解决了收发阵元数大于快拍数时，求解稳健协方差矩阵的过程中产生的秩缺失问题；(3)首次将用于瞬时功率估计中的mw方法应用到大阵列mimo雷达收发角联合估计中，并采用随机矩阵的渐进特征值分布理论作为数学工具，以回波数据样本协方差矩阵的对数形式为基础，基于随机矩阵的围线积分技术和stieltjes变换推导新的代价函数，在收发阵元数和快拍数在同一数量级时具有较好的性能。

本发明提出的方法得到国家自然科学基金项目“基于大维随机矩阵理论的mimo雷达稳健目标检测与估计(no.61371158)”的资助。

附图说明

图1是本发明中所用到的双基地mimo雷达系统模型；

图2是在student-t分布的非高斯噪声背景下基于线性收缩和随机矩阵理论的大阵列双基地mimo雷达目标收发角联合估计方法流程图；

图3是为高斯噪声背景下2d-mw方法、2d-rmw方法和2d-rlsmw三种方法的均方根误差随信噪比变化的曲线，仿真参数设置如下：mn＝36，l＝50；

图4是在student-t分布的非高斯噪声背景下2d-mw方法、2d-rmw方法和2d-rlsmw三种方法的均方根误差随信噪比变化的曲线，仿真参数设置如下：mn＝36，l＝50；

图5(a)是在student-t分布的非高斯噪声背景下利用2d-rmw方法进行收发角联合估计的谱峰搜索图，仿真参数设置如下：mn＝36，l＝30；

图5(b)是图5(a)的等高线图；

图6(a)是在student-t分布的非高斯噪声背景下利用2d-rlsmw方法进行收发角联合估计的谱峰搜索图，仿真参数设置如下：mn＝36，l＝30；

图6(b)是图6(a)的等高线图；

图7在student-t分布的非高斯噪声背景下分别用2d-rmw方法和2d-rlsmw方法进行收发角估计的均方根误差随信噪比变化的仿真曲线，仿真参数设置如下：mn＝36，l＝30。

具体实施方式

图1是双基地mimo雷达的系统模型，m个发射阵元同时发射窄带信号照射k个远场目标，在接收端通过n个接收阵元对回波信号进行综合处理，本发明中发射阵列和接收阵列均采用均匀线阵，阵元间距dt＝dr＝λ/2，发射角度为θ，接收角度为第k个目标的发射角度和接收角度分别为θk和快拍数为l。图2是在student-t分布的非高斯噪声背景下基于线性收缩和随机矩阵理论的大阵列双基地mimo雷达目标收发角估计方法的流程图，其具体实施步骤如下：

步骤一：mimo雷达接收端n个接收阵元对多次快拍得到的回波数据信号进行收集和处理。第l个快拍期间回波数据可以写成一个n×lt的矩阵x^(l)：

其中，表示m行1列的发射方向矢量，t表示发射端；表示n行1列的接收方向矢量，r表示接收端；s表示发射波形矩阵，这里sm(m＝1,…,m)表示第m个发射阵元所发射的波形，其编码长度为lt，且s为归一化正交矩阵，ss^h＝im；表示在第l个快拍期间第k个目标的散射系数；表示一个n行lt列的非高斯噪声矩阵，服从student-t分布；

步骤二：首先对n个接收阵元所接收的信号x^(l)进行匹配滤波处理，得到：

其中，为n行m列的匹配滤波后的信号；

步骤三：对匹配滤波后的结果进行矢量化处理，即y^(l)＝vec(x^(l)s^h)。在第l个快拍期间，经过预处理后的观测数据可表示为

其中，为目标散射系数矩阵；是在第l个快拍期间处理后的噪声矢量，维数为mn行1列；是联合导向矢量矩阵，维数为mn行k列at(θk)和分别为发射和接收方向矢量，是khatri-rao矩阵积，是kronecker矩阵积；

步骤四：收集l个快拍期间的数据y^(l)(l＝1,…,l)，可得到一个观测数据矩阵y，它是一个维数为mn×l的随机矩阵，即：

其中，y＝[y⁽¹⁾,…,y^(l)]，b＝[b⁽¹⁾,…,b^(l)]，n＝[n⁽¹⁾,…,n^(l)]，这里表示包含收发角度的联合导向矩阵，b是服从均值为0，方差为1的目标散射系数矩阵，n为mn×l维的服从student-t分布非高斯噪声矩阵，由高斯噪声矩阵的元素乘上服从伽马分布γ(v,1/v)的因子构成，其中，v＝0.02；

步骤五：接收信号的采样协方差矩阵可以表示为求解协方差矩阵的极大似然估计，对接收信号y^(l)进行归一化yi的概率密度函数为其中u(mn/2)/2π^mn/2为单位圆上的均匀密度分布函数，对p(yi:r)求导并令导数等于零可以得到经过定点迭代可以得到极大似然估计定点迭代的初始矩阵可以是任意正定矩阵，本发明中令为单位矩阵。

步骤六：当收发阵元数大于快拍数时，l<mn，上述极大似然估计方法因为秩缺失而导致无解，此时采用线性收缩技术使得定点迭代的过程能持续进行；收缩过程如下：其中0<ρ<1为收缩系数，将定点迭代完成之后最终收敛的协方差矩阵记为收缩系数ρo由最小均方误差准则得到通过化简可以解得收缩系数为由于r未知，用替代r可以得到

步骤七：在获得稳健协方差矩阵估计的基础上，利用2d-mw方法可以得到收发角联合估计的谱峰搜索函数。mw方法是一种pisarenko结构方法，用于瞬时功率估计，本发明将其应用到双基地mimo雷达目标收发角联合估计中，提出2d-mw方法。pisarenko结构过程如下，对采样协方差矩阵进行特征值分解得到：

其中，(λmn>0)，ee^h＝i，定义0<z<∞上的一个连续的函数f(z)，其反函数为可以得到

收发角联合估计估计的谱峰搜索函数可以表示为

改变函数f(z)可以得到不同的谱峰搜索函数，令函数f(z)＝ln(z)，得到midway(mw)方法的代价函数

步骤八：当收发阵元数与快拍数在同一数量级时，即mn→∞，l→∞，mn/l＝c，c为一个常数，y可以看作一个大维随机矩阵，传统方法性能急剧下降，不能正确估计出收发角度，因此，本发明应用随机矩阵理论中的渐近谱理论和围线积分技术来产生一个稳健的谱峰搜索函数，方法如下：首先根据留数定理，可以表示为

其中，ζ^-是一条在上的顺时针方向曲线，并且该曲线包含r的所有特征值。

对曲线ζ^-进行参数化处理，令其中，是一个包含所有特征值的矩形，表示该矩形的边界，wmn(z)与特征值有关，其关系式如下：

为了简化公式，定义两个变量bmn(z)和mmn(z)来分别描述特征值和特征向量：

得到参数化表示形式为

通过变量bmn(z)，得到wmn(z)的导数对上式进一步化简得到：

当mn→∞，l→∞，mn/l＝c时，用bmn(z)和mmn(z)来描述特征值和特征向量不再准确，此时需要研究特征值和特征向量的渐近表现。在随机矩阵理论中，采样协方差矩阵特征值的瞬时分布函数的stieltjes变换定义如下：

仅仅体现了采样协方差矩阵的特征值渐近特性，为了同时体现特征值和特征向量的渐近特性，定义如下函数：

此时，我们用和替换bmn(z)和mmn(z)，得到：

其中，

为了求解定义一个关于的函数：

其中，x∈[0,c]；

在x∈(0,c)区间上对求导：

利用留数定理可以得到：

其中，是下面等式的解：

可以由柯西留数定理直接得到，再此基础上对求积分，可以得到如下：

令x＝c，经过化简可以得到ξmn(c)：

其中，是下面等式的解：

至此，基于线性收缩和随机矩阵理论的mimo雷达收发角联合估计谱峰搜索函数可以写为：

其中，如上所述；

步骤九：对进行谱峰搜索可以获得发射角θ和接收角

下面通过matlab仿真实验来进一步说明本发明的优点。分别在高斯噪声背景下和student-t分布的非高斯噪声背景下比较了2d-mw、基于随机矩阵理论的2d-rmw、基于线性收缩和随机矩阵理论相结合的2d-rlsmw的性能，验证了本发明所提出方法的有效性。考虑发射阵元数为m＝6、接收阵元数为n＝6，收发阵列为均匀阵列，阵元间距dt＝dr＝λ/2，k＝2个远场目标，目标位置选取快拍数分别取l＝50和l＝30。

图3和图4分别为高斯噪声背景下和student-t分布的非高斯噪声背景下，当收发阵元数和快拍数在同一数量级时，2d-mw方法、2d-rmw方法和2d-rlsmw三种方法的均方根误差随信噪比变化的曲线。仿真条件如下：噪声矩阵n为非高斯噪声，服从student-t分布，矩阵由高斯噪声矩阵的元素乘上服从伽马分布γ(v,1/v)的因子构成，其中，v＝0.02，快拍数l＝50，两个目标的坐标为进行40次蒙特卡罗仿真实验。从图3可以看出，2d-rmw方法和2d-rlsmw方法均方根误差曲线相近，说明在高斯噪声环境下两种方法的估计性能相近；从图4可以看出，相同信噪比下，2d-rlsmw方法的均方根误差要小于2d-rmw方法，说明在student-t分布的非高斯噪声背景下，基于线性收缩和随机矩阵理论的2d-rlsmw方法具有更好的优越性。

图5(a)、图5(b)和图6(a)、图6(b)分别为student-t分布的非高斯噪声背景下收发阵元数大于快拍数时2d-rmw方法和2d-rlsmw方法的谱峰搜索示意图。将快拍数设为l＝30，即收发阵元数大于快拍数的情况，目标位置信噪比为snr＝-5db。从图5(a)和图5(b)可以看出，2d-rmw方法谱峰搜索具有较大的伪峰，无法正确的估计出收发角度；从图6(a)和图6(b)可以看出，2d-rlsmw方法通过谱峰搜索完全可以正确地估计出收发角，充分验证了2d-rlsmw方法的有效性。

图7为相同条件下进行40次蒙特卡罗仿真实验，利用2d-rmw和2d-rlsmw两种方法进行收发角估计的均方根误差(rmse)随信噪比变化的曲线。从图7可以看出，随着信噪比的增加，2d-rmw方法的均方根误差在10°和100°之间浮动，该方法完全失效，无法估计出目标收发角，所估计的收发角度发生随机性变化；2d-rlsmw方法的均方根误差曲线随信噪比的增加呈下降趋势，而且在信噪比为snr＝5db时，角度估计误差小于0.1°，说明在student-t分布的非高斯噪声环境下，当收发阵元数大于快拍数时，基于线性收缩和随机矩阵理论的2d-rlsmw方法具有较大的优势。因此，本发明所提出的方法显著提高了大阵列mimo雷达目标收发角联合估计的稳健性。