一种基于M估计的低快拍数下波达方向估计方法与流程

本发明属于无线移动通信技术领域，具体涉及一种基于m估计的低快拍数下波达方向估计方法。

背景技术：

阵列信号处理经过近几十年的不断发展，已广泛应用于通信、雷达、声呐、电子对抗等许多领域。波达方向估计技术主要研究如何估计空间目标信号的波达方向(directionofarrival，简称doa)，是阵列信号处理领域重要的研究内容之一。无论是在军事还是在民用领域，空间目标信号的波达方向都是极为重要的估计参数，是后续定位、探测、导航和成像等诸多技术的先决基础，因此研究波达方向估计技术有着十分重要的现实意义。

随着信息技术的飞速发展，天线阵列在接收目标发射的信号是受到的干扰日趋多样化和复杂化。在这种背景下，信号接收时是面临的噪声常常不再符合高斯白噪声的特性，而经常呈现为存在“拖尾现象”的复椭圆对称分布。另外，随着现在阵列中天线数目的不断增加，保持大量的快拍数会导致较高的计算复杂度和硬件系统实现的复杂度，因此，如何在小快拍数下保持良好的doa估计精度也是目前阵列信号处理中的一项重要研究。

传统的doa估计算法，如music算法等算法，需要对接收信号的协方差矩阵进行子空间分解来实现。但是在实际应用中，该协方差矩阵是难以直接获得的，因此一般的做法是利用接收信号的样本协方差矩阵代替真实协方差矩阵来进行子空间分解，但是在非高斯噪声的条件下，样本协方差矩阵并不能良好地近似真实的接收信号协方差矩阵，因此基于子空间分解的算法在非高斯噪声的背景下性能会严重降低。最近，一种基于m估计的协方差矩阵估计理论被提出，这种估计方法被证明可以有效对抗存在离群值的“拖尾噪声”。另外，在快拍数不足的情况下，对样本协方差矩阵进行子空间分解得到的特征值和特征向量与真实值偏差很大，因此基于子空间分解的doa估计方法也不能很好地工作。为了解决这个问题，学者们陆续提出了g-music算法、矩阵收缩估计方法和矩阵托普利兹化的方法，但是它们各有缺陷：g-music算法只对空间邻近信号有效、矩阵收缩估计在doa估计中没有合适的收缩目标、矩阵托普利兹化在高信噪比时性能不佳。

另外，在非高斯噪声和快拍数较少两种情况同时存在时，学者们也进行了一定的研究，代表的成果有稳健g-music算法和将矩阵收缩估计与m估计结合的方法，但是前者对空间分离信号没有性能提升，适用范围很窄，后者选取的目标矩阵为单位矩阵，不适用于doa估计领域。因此，本问题目前仍没有得到很好的解决。

技术实现要素：

鉴于以上所述现有技术的缺点，本发明的目的在于提供一种基于m估计的低快拍数下波达方向估计方法，以达到准确实现doa估计的目的。

为实现上述目的及其他相关目的，本发明提供一种基于m估计的低快拍数下波达方向估计方法，该估计方法包括：

步骤1.设置天线阵列并获得所述天线阵列的观测数据矢量x(t)；

步骤2.计算观测数据矢量x(t)在快拍数为t的协方差矩阵

步骤3.对所述样本协方差矩阵进行迹正则化得到迹正则化后的矩阵r1；

步骤4.对迹正则化后的矩阵r1进行托普利兹校正，得到目标矩阵rt；

步骤5.根据迹正则化后的目标矩阵rt，计算矩阵收缩估计的收缩系数ρ；

步骤6.根据所述收缩目标矩阵rt和所述收缩系数ρ，构造新的协方差矩阵r3；

步骤7.将新的协方差矩阵r3重新置为样本协方差矩阵重复步骤2～步骤6，直到得到最终的估计；

步骤8.利用得到的协方差矩阵r3采用子空间分解的方法计算信号的波达方向。

可选地，所述观测数据矢量x(t)的计算方法为：

x(t)＝as(t)+n(t)

其中，x(t)＝[x0(t),…,xm-1(t)]^t为观测数据矢量，s(t)＝[s1(t),…,sk(t)]^t为信号源数据矢量，a＝[a(θ1),…,a(θk)]为阵列流型矩阵，a(θk)＝[a0,k,…,am-1,k]^t为导向矢量，n(t)＝[n0(t),…,nm-1(t)]^t为噪声矢量。

可选地，所述样本协方差矩阵的计算公式为：

其中，x(tn)表示在第tn个快拍时的阵列观测数据矢量。

可选地，所述迹正则化后的矩阵r1通过以下方法计算得到：

其中，trace表示计算矩阵的迹，即矩阵对角线元素的和。

可选地，所述目标矩阵rt通过以下方法计算得到：

其中，j^m代表m维的移位矩阵，j^m中只有第m条对角线元素为1，其余元素为0。

可选地，所述收缩系数ρ通过以下方法计算得到：

可选地，所述新的协方差矩阵r3通过以下方法计算得到：

其中，xi代表观测数据矢量x(t)的第i列，其代表阵列在第i次快拍中的全部观测数据。

可选地，所述步骤8包括以下子步骤：

对新构造的协方差矩阵r3进行特征分解，得到其噪声子空间un和信号子空间us，其中，σs表示对协方差矩阵进行特征值分解后，由较大的k个特征值构成的对角矩阵，σn表示由较小的m-k个特征值成的对角矩阵；

由music算法，利用所述信号子空间us和所述噪声子空间un的正交性得到空间谱函数为：

其中，θ为谱峰搜索时的角度，a(θ)代表阵列的导向矢量，un代表对协方差矩阵进行特征值分解后，由较小的m-k个特征值对应的特征向量构成的矩阵，

如上所述，本发明的一种基于m估计的低快拍数下波达方向估计方法，具有以下有益效果：

本发明使用托普利兹校正的方法克服了快拍数不足导致doa估计精度下降的问题；在构造协方差矩阵时使用了m估计的方法减轻了非高斯噪声对doa估计精度的影响；同时将托普利兹校正后的矩阵作为目标矩阵，m估计得到的矩阵作为原始矩阵，利用协方差矩阵收缩估计的方法构造新的协方差矩阵，克服了托普利兹校正带来的在doa估计在高信噪比下性能不佳的问题；综上，本发明可以在小快拍数和非高斯噪声同时存在的条件下实现准确的doa估计。

附图说明

为了进一步阐述本发明所描述的内容，下面结合附图对本发明的具体实施方式作进一步详细的说明。应当理解，这些附图仅作为典型示例，而不应看作是对本发明的范围的限定。

图1为本发明阵列设置示意图；

图2为本发明实施方式仿真实验doa的均方根误差在高斯白噪声下随snr变化关系示意图；

图3为本发明实施方式仿真实验doa的均方根误差在student-t分布的噪声下随snr变化关系示意图；

图4为本发明实施方式仿真实验doa的均方根误差在威布尔分布的噪声下随snr变化关系示意图；

图5为本发明实施方式仿真实验doa的均方根误差在k分布的噪声下随snr变化关系示意图；

图6为本发明实施例的一种基于m估计的低快拍数下波达方向估计方法的流程图。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，遂图式中仅显示与本发明中有关的组件而非按照实际实施时的组件数目、形状及尺寸绘制，其实际实施时各组件的型态、数量及比例可为一种随意的改变，且其组件布局型态也可能更为复杂。

如图6所示，本发明提供一种基于m估计的低快拍数下波达方向估计方法，包括以下步骤：

步骤1：设置天线阵列并获得所述天线阵列的观测数据矢量x(t)。

具体地，设置一个如图1所示的一维均匀线阵，共有m个阵元，阵列间距d＝λ/2。该阵列接收k个波长为λ的窄带平稳独立且非零峰度信号，k≤m-1，阵列采样的快拍数t略大于阵元数m。在本实施方式中阵元数m设置为20，阵列的快拍数t设置为25。本实施方式有k＝3两个信号源，其入射角度分别为θ1＝-3°，θ2＝16°，θ3＝55°。阵列阵元上的噪声为复椭圆对称分布的非高斯噪声，且噪声独立于信号。

第m个阵元在t时刻接收到的数据可以表示为：

其中，sk(t)表示第k个信号源，nm(t)表示第m个阵元上的噪声，而am,k表示第m个阵元对第k个信号的阵列流型元素，可以进一步表示为：

其中，ωk有下列形式：

其中，λ是信号波长，θk是第k个信号的doa。

将阵列观测数据表示成矩阵形式为：

x(t)＝as(t)+n(t)

其中，x(t)＝[x0(t),…,xm-1(t)]^t为观测数据矢量，s(t)＝[s1(t),…,sk(t)]^t为信号源数据矢量，a＝[a(θ1),…,a(θk)]为阵列流型矩阵，a(θk)＝[a0,k,…,am-1,k]^t为导向矢量，n(t)＝[n0(t),…,nm-1(t)]^t为噪声矢量。

步骤2：设置天线阵列并获得所述天线阵列的观测数据矢量x(t)；

在实际工作中，观测数据矢量x(t)的样本协方差矩阵的计算公式为：

其中，t为快拍数，表示在第tn个快拍时的阵列观测数据矢量。

步骤3：利用样本协方差矩阵计算迹正则化后的矩阵r1，计算公式为：

其中，trace代表计算矩阵的迹，即矩阵对角线元素的和。

步骤4：对迹正则化后的矩阵r1进行托普利兹校正，得到目标矩阵rt；

阵列接收数据的真实协方差矩阵是托普利兹矩阵，但当实际接收信号的快拍数较少时，观测数据的样本协方差矩阵不再满足矩阵的托普利兹特性，因此对步骤3中得到的迹正则化后的矩阵r1进行托普利兹校正，得到目标矩阵rt，以减轻低快拍数带来的影响，计算公式为：

其中，j^m代表m维的移位矩阵，j^m中只有第m条对角线元素为1，其余元素为0。

步骤5：根据迹正则化后的矩阵，计算矩阵收缩估计的收缩系数ρ；

由于托普利兹化使阵列协方差矩阵的真实元素发生了变化，因此信噪比较高时，利用托普利兹化后的矩阵计算的doa的精度会比直接利用样本协方差矩阵计算的doa的精度低。为了解决这个问题，本发明采用矩阵收缩估计的方法，将托普利兹校正后的目标矩阵rt作为收缩目标矩阵，以提升在高信噪比下的性能，采用收缩估计的方法需要计算矩阵收缩系数，本发明采用“障碍估计”的原则计算收缩系数ρ，其计算公式为：

步骤6：根据所述收缩目标矩阵rt和所述收缩系数ρ，构造新的协方差矩阵r3；

本发明采用m估计的思想重构样本协方差矩阵以对抗非高斯噪声，当阵列的观测矩阵x(t)已知时，m估计的表达式为：

其中xi代表观测数据x(t)的第i列，其代表阵列在第i次快拍中的全部观测数据。u要求是非负、非增、连续的函数。本发明采用tyler提出的估计法，即取u(x)＝1/x。另外，在表达式的左右两边都出现了协方差矩阵r，因此该方法需要一个初始的协方差矩阵r作为初值，然后采取迭代的方式进行求解。另外，为了同时对抗小快拍数带来的影响，本发明同时采用协方差矩阵收缩估计的思路，将m估计得到的协方差矩阵作为收缩估计的主矩阵，将步骤4中托普利兹化后得到的矩阵rt作为收缩目标矩阵，构造新的协方差矩阵，其计算公式为：

其中，收缩系数ρ由步骤5中的计算得到。

步骤7：将新的协方差矩阵r3重新置为样本协方差矩阵重复步骤2～步骤6，直到得到最终的估计；

如步骤6所述，m估计是一个迭代的过程，因此在构造出新的协方差矩阵后，我们仍需要对新构造的协方差矩阵r3进行迭代更新以取得更精确的协方差矩阵估计结果。具体而言，将步骤5中得到的协方差矩阵r3重新置为样本协方差重复步骤2～步骤6，直到得到最终的估计r3，在本发明中，迭代次数设置为100。

步骤8：利用得到的协方差矩阵r3采用子空间分解的方法计算信号的波达方向。

对步骤7所得新构造的协方差矩阵r3进行特征分解，得到其噪声子空间un和信号子空间us。并利用最大的3(k)个特征值所对应的特征向量uz1,uz2,uz3作为列建立特征向量矩阵uzs＝{uz1,uz2,uz3}作为信号子空间，利用其余17(m-k)个小特征值对应的特征向量uz4,...,uz20为列建立特征向量矩阵uzn＝{uz4,...,uz20}作为噪声子空间，3个大特征值组成对角矩阵ds，其余17个小特征值组成对角矩阵dn，如下式所示：

再利用信号子空间和噪声子空间的正交性得到music算法的空间谱函数为：

其中，θ为谱峰搜索时的角度。

设置第一次谱峰搜索的步长为l1＝1，由music谱函数得到3个峰值，设为θ'1，θ'2，θ'3；设置第二次搜索范围为θ'1±2，θ'2±2，θ'3±2，搜索步长为l2＝0.1，得到最终的doa的估计值

为验证算法的性能，本实施方式设计四组仿真实验。仿真选取了四种典型的噪声分布场景，一种是高斯白噪声，其余3种均属于在阵列信号处理和雷达领域常见的非高斯噪声，分布为student-t分布的噪声、威布尔分布的噪声和k分布的噪声。其中student-t分布的自由度设置为2.5；威布尔分布的尺度参数设置为1，形状参数设置为0.5；k分布的尺度参数设置为2，形状参数设置为0.25。在仿真实验中，将建议方法与5种其它的方法的性能进行了对比，其余的方法分别为：经典music算法、仅使用了m估计的方法、稳健g-music算法、同时结合托普利兹化与m估计的方法和将m估计与矩阵收缩估计结合的方法。

两组实验的随机实验次数均为500，所用均匀线阵相同，阵元数为20，阵列快拍数为25，阵列间距为d＝λ/2。入射的信号源数k＝3。信号参数分别为θ1＝-3°，θ2＝16°，θ3＝55°。四组实验的结果分别如图2、3、4、5所示。由图3、4、5可知，在非高斯分布的噪声存在的条件下，本发明所提方案在信噪比较低的条件下比其它的方法估计精度更高，同时在信噪比较高的条件下不会出现托普利兹化带来的精度损失问题。

上述实施例仅例示性说明本发明的原理及其功效，而非用于限制本发明。任何熟悉此技术的人士皆可在不违背本发明的精神及范畴下，对上述实施例进行修饰或改变。因此，举凡所属技术领域中具有通常知识者在未脱离本发明所揭示的精神与技术思想下所完成的一切等效修饰或改变，仍应由本发明的权利要求所涵盖。