基于约束张量分解的mimo雷达二维波达方向估计方法
技术领域
本发明属于雷达
技术领域
，具体涉及基于约束张量分解的MIMO雷达二维波达方向估计方法。背景技术
MIMO雷达是一种发射波形相互正交的新体制雷达，因其在多目标检测和参数估计等方面的优越性能，在过去二十年里受到了广泛的关注，其中尤以天线单元间距相对较小的集中式MIMO雷达为代表。众多学者研究并分析了MIMO雷达领域较之传统的相控阵雷达的优势和性能改善，这包括：良好的抗干扰性能，灵活的发射方向图以及分辨率更好、分辨精度更高的波达方向估计。这种性能改善主要来自对波形分集的有效利用，即通过接收端的匹配滤波得到全部收发通道的回波从而等效构建一个口径更大阵元更多的虚拟阵列。正因如此，过去对MIMO雷达波达方向估计的研究主要集中在对虚拟阵列接收数据的协方差矩阵进行分析，这类方法可以看作是对相控阵雷达波达方向估计方法的推广，比如[Z.Guo,X.Wang,andW.Heng,“Millimeter-wave channel estimation basedon2-D beamspaceMUSIC method,”IEEE Trans.Wireless Commun.,vol.16,no.8,pp.5384–5394,2017]，其将传统的MUSIC算法推广到MIMO雷达；[C.Jinli,G.Hong,and S.Weimin,“Angle estimationusing ESPRIT withoutpairing in MIMO radar,”Electron.Lett.,vol.44,no.24,pp.1422–1423,2008.]，则是考虑了ESPRIT算法在MIMO雷达的应用。这些算法一次仅能利用MIMO雷达多脉冲接收数据中的单一脉冲数据，并且需要在不同脉冲之间迭代波达方向估计结果，既容易受到目标波动的影响，在目标回波信噪比较低的时候又不能保持较好的估计效果。
针对以上问题，有学者提出基于张量分解的MIMO雷达波达方向估计算法，比如[D.Nion and N.D.Sidiropoulos,“Tensor algebra and multidimensional harmonicretrieval in signal processing for MIMO radar,”IEEE Trans.Signal Process.,vol.58,no.11,pp.5693–5705,Nov.2010.]和[N.D.Sidiropoulos,L.De Lathauwer etal.,“Tensor decomposition for signal processing and machine learning,”IEEETrans.Signal Process.,vol.65,no.13,pp.3551–3582,Jul.2017.]。采用张量模型同时存储MIMO雷达多脉冲接收数据，既可以利用MIMO雷达接收数据之间的多线性结构，又能够同时对多个目标的波达方向进行估计，且有效改善了波达方向估计的性能。然而，常规的张量分解方法，即交替最小二乘法(Alternating Least Squares)，计算复杂度高、收敛性不稳定、要求目标的数量信息作为先验条件。这些问题在目标张量高于三阶的时候更加明显。
在某些应用场合，比如针对多目标检测跟踪的地基雷达，发射阵列一般具有大量阵元，为了简化系统结构，会采取子阵划分等手段，其对应的张量模型可能达到四阶甚至更高。此时，要对多个目标进行近似实时的波达方向估计，要求所采用的波达估计算法：
因此，目前亟需一种雷达估计算法，能够适用于高阶张量模型，在目标数量未知情况下仍然能够有效进行估计，计算复杂度低，收敛性稳定，波达方向估计精度高且分辨率好。
发明内容
有鉴于此，本发明提供基于约束张量分解的MIMO雷达二维波达方向估计方法，该方法计算复杂度低，收敛性稳定，在目标数未知的情况下仍然能够保持有效估计，可以应用于具备多个相同子阵的MIMO雷达场景，且其角度估计的分辨率和准确度较之其他算法有显著改善，为MIMO雷达工程应用提供技术途径。
为达到上述目的，本发明的技术方案为：
基于约束张量分解的MIMO雷达二维波达方向估计方法，包括如下步骤：
(1a)构造用于MIMO雷达接收数据四阶张量模型，MIMO雷达具有多个相同发射子阵。
(1b)重构四阶张量模型并对获得的重构矩阵做奇异值分解，分解结果包含所求的左奇异矩阵，左奇异矩阵包含范德蒙德结构。
(1c)根据左奇异矩阵的范德蒙德结构，利用左奇异矩阵的子矩阵的相互关系，估计特征值的列向量并根据该列向量进行计算，得到二维波达方向的估计值。
进一步的，构造用于MIMO雷达接收数据四阶张量模型，MIMO雷达具有多个相同发射子阵，具体方法为：
(2a)构建MIMO雷达的接收阵列和第s个发射子阵As在第q个脉冲的接收数据其中s＝1，2,...S，q＝1，2,...，Q，S＝I×J，S表示总发射子阵数目，I表示发射阵列在x轴方向的子阵数目，J表示发射阵列在y轴方向的子阵数目，Q表示总脉冲数，每个发射子阵有个阵元均匀分布在矩阵网格上，的表达式为：
其中，为M0的子阵列在x轴方向的发射阵元数目，为M0的子阵列在x轴方向的发射阵元数目，M0为每个发射子阵所包含的阵元数目，表示接收阵列和第s个发射子阵在第q个脉冲的回波数据，As表示第s个发射子阵的发射导向矩阵，B表示接收阵列的导向矩阵，∑q＝diag(cq)是由列向量cq张开的方阵，包含L个目标的多普勒频移和雷达反射系数信息， 代表目标多普勒矢量，σ反应目标雷达反射系数，对l＝1,2,...L，L为总目标数，fl表示第l个目标的多普勒频移，T是雷达的脉冲重复周期，是对应的高斯白噪声矩阵，(.)T表示矩阵转置，表示定义为的意思。
(2b)利用发射阵列均匀排布结构，分析各发射子阵的发射导向矢量之间的关系，第s个发射子阵As、As对应的横向导向矩阵分量Uj和纵向导向矩阵分量Vi的表达式为：
AS＝Uj⊙Vi
Uj＝U0Γj,Vi＝V0Γi
其中，s＝(j-1)I+i表示第s个发射子阵的序号，j＝1,2,...J和i＝1，2，...I分别表示发射子阵在横向和纵向的编号，Uj是发射子阵的横向导向矩阵分量，Vi是发射子阵的纵向导向矩阵分量，U0是横向参考矩阵，V0是纵向参考子阵，Γj是包含i个列向量张开的对角阵，Γi是由包含j个列向量张开的对角阵；⊙表示Khatri-Rao积。
(2c)将列向量化，得到 A0＝U0⊙V0是参考发射子阵的导向矩阵，将高斯白噪声矩阵列向量化，得到并将步骤(2b)中描述的关系式代入其中，按照子阵编号顺序将全部S个列向量组合成一个新的矩阵Y(q)，即表达为：
其中，代表各发射子阵在横向的相位信息集合， 代表各发射子阵在纵向的相位信息集合， 是阵列在第q个脉冲接收到的噪声矩阵。
(2d)列向量化矩阵Y(q)得到zq，具体写作zq＝[H⊙Δ⊙A0⊙B]cq+rq，rq是N(q)的列向量化结果，然后按列拼接Q个脉冲的接收数据得到MIMO雷达多脉冲接收数据Z，即表达式为：
Z＝[H⊙Δ⊙A0⊙B]CT+R
其中是Q个脉冲的目标多普勒和RCS信息集合， 表示阵列在Q个脉冲接收的全部噪声集合。
(2e)将矩阵Z重构为一个5阶张量表达式为：
其中，表示向量外积，[[.]]用来表示张量矩阵因子的集合，ηl，δl，αl，βl，γl分别是H，Δ，A0，B，C的第l列向量，是由矩阵R经过相同重构方式生成的5阶噪声张量。
进一步的，重构四阶张量模型并对获得的重构矩阵做奇异值分解，分解结果包含所求的左奇异矩阵，左奇异矩阵包含范德蒙德结构，具体方法为：
(3a)将张量通过张量变换化为一个三阶张量表达式为：
其中，G＝H⊙Δ是第一矩阵因子，第三矩阵因子B⊙C满足列满秩条件。
(3b)根据张量重构原则，将三阶张量从第三个维度进行重构，得到重构后的矩阵T(3),表达式为：
T(3)＝(G⊙A0)(B⊙C)T+N
其中，矩阵N表示由噪声张量通过相同重构方式产生的噪声矩阵。
(3c)对重构后的矩阵T(3)做奇异值分解，表达式为：
T(3)＝U∑VH
其中，(.)H表示矩阵共轭转置，分解结果分别是维度大小为SM0×L的左奇异矩阵U，维度大小为NQ×L的右奇异矩阵V，和维度为L×L的奇异值矩阵∑；由于第三因子矩阵B⊙C列满秩，必然存在一个L×L的非奇异变换矩阵E，使得UE＝H⊙Δ⊙A0，而矩阵H和Δ都是范德蒙德矩阵。
进一步的，根据左奇异矩阵的范德蒙德结构，利用左奇异矩阵的子矩阵的相互关系，估计特征值的列向量并根据该列向量进行计算，得到二维波达方向的估计值，具体方法为：
(4a)定义左奇异矩阵的两个子矩阵第一子矩阵U1和第二子矩阵U2，表达式为：
其中，是由矩阵H除第一行外组成的子矩阵，是由矩阵H除最后一行外组成的子矩阵，矩阵E为非奇异变换矩阵。
(4b)根据Khatri-Rao积的运算规律，利用行选择从矩阵U中构造子矩阵U1和U2，表达式为：
其中，I表示单位矩阵，其维度大小由对应下标决定，0表示一个全部元素为0的矩阵，其维度大小也由对应下标决定。
(4c)利用范德蒙德矩阵结构的特点，U1和U2的关系表达式为：
U2E＝U1EΩy
其中，Ωy＝diag(ωy)是由ωy张开的方阵，是范德蒙德矩阵H对应的生成因子向量，u1为ul的第1个导向矢量，ul为矩阵U的导向矢量的第一中间量，Δmy＝mj+1-mj表示第j+1个发射子阵和第j个发射子阵在横向的相位步进。
(4d)定义左奇异矩阵的另外两个子矩阵：第三子矩阵U3和第四子矩阵U4，表达式为：
其中，是由矩阵Δ除第一行外组成的子矩阵，是由矩阵Δ除最后一行外组成的子矩阵。
(4e)根据Khatri-Rao积的运算规律，利用行选择从矩阵U中构造第三子矩阵U3和第四子矩阵U4，表达式为：
其中，表示Kronecker积。
(4f)利用范德蒙德矩阵结构的特点，U3和U4的关系表达式为：
U4E＝U3EΩx
其中，Ωx＝diag(ωx)是由ωx张开的方阵，是范德蒙德矩阵Δ对应的生成因子向量，Δmx＝mi+1-mi第i+1个发射子阵和第i个发射子阵在纵向的相位步进。
(4g)从各子矩阵的相互关系中，得到对范德蒙德矩阵H对应的生成因子向量ωy和范德蒙德矩阵Δ对应的生成因子向量ωx的估计，表达式为：
其中，(.)-1表示矩阵求逆；分别对矩阵和做特征分解，得到的特征值所组成的列向量记作和即为对列向量ωy和ωx的估计。
(4h)根据和的值，实现目标的二维波达方向估计，表达式为：
其中，和分别是列向量和的第l个元素，和表示对中间量的估计，而和即为对第l个目标俯仰角和方位角的二维波达方向估计结果；收集全部L组估计值即实现所提二维波达方向估计。
有益效果：本发明方法根据MIMO雷达多脉冲接收数据的多重线性结构构造高阶张量模型及其对应的矩阵展开，对展开矩阵做奇异值分解并通过因发射阵列结构带来的约束条件，从左奇异矩阵中产生四个满足旋转不变性的子矩阵，随后利用一种类ESPRIT算法对多个目标进行二维波达方向估计。本发明方法利用张量模型，挖掘了MIMO雷达多脉冲接收数据之间的多重线性关系；不依赖目标数目的先验信息，能够在其未知的情况下有效估计出目标个数，并同时对多目标进行波达方向估计；仅涉及矩阵运算，不要求迭代，因此计算复杂度低且一定保证收敛，无噪声情况下其分解结果与原始输入条件完全一致。同时，本发明方法对波达方向估计的分辨率更高、分辨精度更好。
附图说明
图1是本发明的系统阵列结构示意图；
图2是本发明的算法流程图；
图3是本发明仿真合成的波达方向估计精度示意图；
图4是本发明仿真合成的波达方向估计分辨率示意图。
具体实施方式
下面结合附图并举实施例，对本发明进行详细描述。
本发明的核心思路是：首先推导具有多个相同发射子阵的MIMO雷达接收数据的高阶张量模型，并论证了其对应的因子矩阵受阵列结构约束具有范德蒙德(Vandermonde)结构,然后通过对该高阶张量的矩阵重构和奇异值分解，保证左奇异矩阵继承了范德蒙德结构，最后利用一种类ESPRIT算法，从左奇异矩阵的子矩阵中估计出目标来向的相位信息，进而实现二维波达方向估计。
本发明的技术原理是：如图1所示，考虑MIMO雷达具有M＝MxMy个均匀分布在二维矩形网格上的发射阵元，其中Mx和My分别表示该矩形网格在x轴和y轴方向上的点数。相邻阵元之间的间距为其中λ是雷达系统的工作波长。不失一般性，假设电磁波在空间中传播无衰减且满足远场条件，则该阵列在俯仰角方位角θ方向的发射导向矢量可以表示为
其中， (·)T表示矩阵转置，表示定义一个变量。类似地，假定该雷达具有N个从发射阵元中随机选取得到的接收阵元，其坐标为(xn,yn),n＝1,2,...,N，则相应的接收导向矢量可以表示为
一般而言，当空间中存在L个目标，其分布为和可以分别定义发射导向矩阵和接收导向矩阵为
如此，该MIMO雷达的多目标接收数据可以表示为
Y＝BΣAT+N
其中，Σ＝diag(σ)，表示L个目标雷达反射系数(RCS)的集合，N表示一个N×M的高斯白噪声矩阵。二维波达方向估计，就是从对Y的观测中获取全部目标的角度信息。如图3所示，该图是本发明仿真合成的波达方向估计精度示意图，图4是本发明仿真合成的波达方向估计分辨率示意图。由两图可知，本发明方法能够提高波达方向的估计定都和估计分辨率。
如图2所示，本发明提出一种基于约束张量分解的MIMO雷达二维波达方向估计方法，包括如下步骤：
(1a)对于一个具有M个分布在矩形网格的发射阵元和具有N个随机分布的接收阵元的MIMO雷达，构造用于MIMO雷达接收数据四阶张量模型，该MIMO雷达具有多个相同的发射子阵：
(2a)构建MIMO雷达的接收阵列和第s个发射子阵As在第q个脉冲的接收数据其中s＝1，2,...S，q＝1，2,...，Q，S＝I×J，S表示总发射子阵数目，I表示发射阵列在x轴方向的子阵数目，J表示发射阵列在y轴方向的子阵数目，Q表示总脉冲数。其中，矩形网格的长为x轴，矩形网格的宽为y轴，矩形网格的第一个元素为原点。每个发射子阵有个阵元均匀分布在矩阵网格上，的表达式为
其中，为M0的子阵列在x轴方向的发射阵元数目，为M0的子阵列在x轴方向的发射阵元数目，M0为每个发射子阵所包含的阵元数目，表示接收阵列和第s个发射子阵在第q个脉冲的回波数据，As表示第s个发射子阵的发射导向矩阵，B表示接收阵列的导向矩阵，Σq＝diag(cq)是由列向量cq张开的方阵，包含L个目标的多普勒频移和雷达反射系数信息， 代表目标多普勒矢量，σ反应目标雷达反射系数，对l＝1，2，...L，L为总目标数，fl表示第l个目标的多普勒频移，T是雷达的脉冲重复周期，是对应的高斯白噪声矩阵，(.)T表示矩阵转置，表示定义为的意思。
(2b)利用发射阵列均匀排布结构，分析各发射子阵的发射导向矢量之间的关系，第s个发射子阵As、As对应的横向导向矩阵分量Uj和纵向导向矩阵分量Vi的表达式为：
As＝Uj⊙Vi
Uj＝U0Γj，Vi＝V0Γi
其中，s＝(j-1)I+i表示第s个发射子阵的序号，j＝1，2，...J和i＝1，2，...I分别表示发射子阵在横向和纵向的编号，Uj是发射子阵的横向导向矩阵分量，Vi是发射子阵的纵向导向矩阵分量，U0是横向参考矩阵，V0是纵向参考子阵，Γj是包含i个元素的列向量张开的对角阵，Γi是由包含j个元素的列向量张开的对角阵。本发明中，选取第一个发射子阵作为参考矩阵。⊙表示Khatri-Rao积。
此外，U、V的表达式为：
收集全部L组因发射阵列产生的导向矢量，右奇异矩阵V的导向矢量写作：
左奇异矩阵U的导向矢量写作：
第一中间量为矩阵U的第l个导向矢量，第二中间量 为矩阵V的第l个导向矢量。两者包含目标的方位角信息θl和俯仰角信息Γj＝diag(hj)为列向量hj张开的对角阵，Γi＝diag(di)为列向量di张开的对角阵，该组列向量代表子阵间的相位信息，对第s个子阵相位中心坐标(mi，mj)，两个对角阵的列向量具体表示为
(2c)将列向量化，得到 A0＝U0⊙V0是参考发射子阵的导向矩阵，将高斯白噪声矩阵列向量化，得到并将步骤(2b)中描述的关系式代入其中，按照子阵编号顺序将全部S个列向量组合成一个新的矩阵Y(q)，即表达为：
其中，代表各发射子阵在横向的相位信息集合， 代表各发射子阵在纵向的相位信息集合， 是阵列在第q个脉冲接收到的噪声矩阵。
(2d)列向量化矩阵Y(q)得到zq，具体写作zq＝[H⊙Δ⊙A0⊙B]cq+rq，rq是N(q)的列向量化结果，然后按列拼接Q个脉冲的接收数据得到MIMO雷达多脉冲接收数据Z，即所述表达式为：
Z＝[H⊙Δ⊙A0⊙B]CT+R
其中是Q个脉冲的目标多普勒和RCS信息集合， 表示阵列在Q个脉冲接收的全部噪声集合。
(2e)将矩阵Z重构为一个5阶张量所述表达式为：
其中，表示向量外积，[[.]]用来表示张量因子矩阵的集合，ηl，δl，αl，βl，γl分别是H，Δ，A0，B，C的第l列向量，是由矩阵R经过相同重构方式生成的5阶噪声张量。
(1b)重构所述四阶张量模型并对获得的重构矩阵做奇异值分解，分解结果包含所求的左奇异矩阵，该左奇异矩阵包含范德蒙德结构：
(3a)将张量通过张量变换化为一个三阶张量表达式为：
其中，G＝H⊙Δ是第一矩阵因子，A0可作为第二矩阵因子使用，第三因子矩阵B⊙C满足列满秩条件。
(3b)根据张量重构原则，将三阶张量从第三个维度进行重构，得到重构后的矩阵T(3),表达式为：
T(3)＝(G⊙A0)(B⊙C)T+N
其中，矩阵N表示由噪声张量通过相同重构方式产生的噪声矩阵。
(3c)对重构后的矩阵T(3)做奇异值分解，表达式为：
T(3)＝U∑VH
其中，(.)H表示矩阵共轭转置，分解结果分别是维度大小为SM0×L的左奇异矩阵U，维度大小为NQ×L的右奇异矩阵V，和维度为L×L的奇异值矩阵∑。由于第三因子矩阵B⊙C列满秩，必然存在一个L×L的非奇异变换矩阵E，使得UE＝H⊙Δ⊙A0，而矩阵H和Δ都是范德蒙德矩阵。上式表明，左奇异矩阵受到阵列结构的约束存在特殊的范德蒙德结构。
(1c)根据左奇异矩阵的范德蒙德结构，利用左奇异矩阵的子矩阵的相互关系，估计特征值的列向量并根据该列向量进行计算，得到二维波达方向的估计值。
(4a)定义左奇异矩阵的两个子矩阵第一子矩阵U1和第二子矩阵U2，表达式为：
其中，是由矩阵H除第一行外组成的子矩阵，是由矩阵H除最后一行外组成的子矩阵，矩阵E为非奇异变换矩阵。
(4b)根据Khatri-Rao积的运算规律，利用行选择从矩阵U中构造子矩阵U1和U2，表达式为：
其中，I表示单位矩阵，其维度大小由对应下标决定，0表示一个全部元素为0的矩阵，其维度大小也由对应下标决定。
(4c)利用范德蒙德矩阵结构的特点，U1和U2的关系表达式为：
U2E＝U1EΩy
其中，Ωy＝diag(ωy)是由ωy张开的方阵，是范德蒙德矩阵H对应的生成因子向量，u1为ul的第1个导向矢量，ul为矩阵U的导向矢量的第一中间量，Δmy＝mj+1-mj表示第j+1个发射子阵和第j个发射子阵在横向的相位步进。
(4d)定义左奇异矩阵的另外两个子矩阵U3和U4，表达式为：
其中，是由矩阵Δ除第一行外组成的子矩阵，是由矩阵Δ除最后一行外组成的子矩阵。
(4e)根据Khatri-Rao积的运算规律，利用行选择从矩阵U中构造子矩阵U3和U4，表达式为：
其中，表示Kronecker积。
(4f)利用范德蒙德矩阵结构的特点，U3和U4的关系表达式为：
U4E＝U3EΩx
其中，Ωx＝diag(ωx)是由ωx张开的方阵，是范德蒙德矩阵Δ对应的生成因子向量，表示第i+1个发射子阵和第i个发射子阵在纵向的相位步进。
(4g)从各子矩阵的相互关系中，得到对范德蒙德矩阵H对应的生成因子向量ωy和范德蒙德矩阵Δ对应的生成因子向量ωx的估计，表达式为：
其中，(.)-1表示矩阵求逆；分别对矩阵和做特征分解，得到的特征值所组成的列向量记作和即为对列向量ωy和ωx的估计。
(4h)根据和的值，实现目标的二维波达方向估计，表达式为：
其中，和分别是列向量和的第l个元素，和表示对中间量的估计，而和即为对第l个目标俯仰角和方位角的二维波达方向估计结果。收集全部L组估计值即可实现所提二维波达方向估计。
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
- 上一篇：石墨接头机器人自动装卡簧、装栓机
- 下一篇：一种多雷达航迹融合方法
