基于张量稀疏表示的MIMO雷达成像方法与流程

本发明属于雷达技术领域和信号处理领域，特别涉及一种多输入多输出类型系统的基于张量稀疏表示的MIMO雷达成像方法。

背景技术：

多输入多输出(MIMO)雷达是一种21世纪新出现的雷达系统，它利用多个发射与接收天线同时对目标进行观测。良好的阵列构型设计和波形分集技术使得MIMO雷达能够获得远多于实际物理阵元个数的观测通道和空间自由度，可以显著的改善参数的可辨识性，实现更为灵活的发射方向图设计，改进目标检测和参数估计性能。相比于传统成像雷达，MIMO雷达在成像的方位向分辨率、实时性和运动补偿方面具有明显的性能优势。因此MIMO雷达成像具有广泛的应用前景。

常见的MIMO雷达成像方法，诸如BP(back projection)方法或DAS(delay and sum)类波束形成方法，包括改进的Kirchhoff偏移方法、衍射堆栈方法等，具有与匹配滤波和波束形成相似的形式，其优点是方法简单易于实现，输出信噪比高，但是存在分辨率较低且旁瓣水平高，成像效果差的缺陷。

为了获得更好的成像效果，人们将压缩感知技术应用到MIMO雷达成像中。稀疏微波成像是指将压缩感知与雷达成像有机结合形成的一种新的成像方法。它通过寻找被观测目标的少量回波数据，利用稀疏重构技术提取目标的空间位置、散射特征和运动特征等参数。和传统的雷达成像方法相比，压缩感知的引入可以显著地降低系统的数据采集率和系统复杂度，而且稀疏重构方法潜在的超分辨能力有进一步提升成像性能的能力。在压缩感知框架下，MIMO雷达成像可以视为一个稀疏估计问题，其成像过程可以用线性规划的方法或是贪婪类方法予以解决。佛罗里达大学的Li J.教授等提出了许多适用于MIMO雷达稀疏成像的稀疏重构方法，如循环自适应方法(Iterative Adaptive Approaches,IAA)和稀疏学习循环最小化方法(Sparse Learning via Iterative Minimization,SLIM)等。Higgins等提出了空间-距离自适应处理方法。这些MIMO雷达成像方法都是将自适应技术应用到二维联合滤波器权矢量的设计中，通过迭代更新二维权矢量和获得的图像幅值，通过一定的迭代次数最终得到高分辨率和低旁瓣的成像结果。但是，这些方法自适应维数巨大，方法的时间复杂度过高，不仅难以进行实时成像，而且在常用的处理器上运行都极为困难。文献(Joint wall mitigation and compressive sensing for indoor image reconstructon.IEEE Transaction on Geosci.Remote Sensing,2013,51(2):891-906.)采用线性规划的方法解决压缩感知问题，得到良好的效果。贪婪恢复方法的代表是OMP类方法。这种方法有着较低的运算负荷，较高的成像分辨率，但是由于OMP方法在基信号选择时只能扩充不能去除不良基信号的策略，OMP类恢复方法在雷达成像应用中会存在伪影点，这不利于目标的识别。文献(Subspace pursuit for compressive sensing signal recomstruction.IEEE Transaction on Information Theory,2009,55(5):2230-2249)提出了被称为子空间追踪方法(subspace pursuit,SP)的压缩感知贪婪方法，纠正了OMP方法中存在伪影点的问题，但是在MIMO雷达成像应用中其分辨率较OMP方法低。

此外，上述所提的稀疏成像方法通常将接收信号作堆栈处理，这无疑将破坏信号的多维结构从而不能利用信号的多维结构信息导致方法性能降低。特别在低信噪比和小样本情况下，性能更加恶化。本发明提出将张量信号处理方法引入MIMO雷达稀疏成像中去。

技术实现要素：

本发明的目的在于提供一种基于张量稀疏表示的MIMO雷达成像方法。

本发明的目的是这样实现的：

本发明包括如下步骤：

(1)M个发射阵元发射相互正交的相位编码信号，N个接收阵元接收该相位编码信号；

(2)应用匹配滤波器对接收到的雷达信号进行匹配滤波；

(3)对匹配滤波后的信号做傅里叶变换，得到空间谱域回波表达式；

(4)场景进行网格划分，将雷达回波离散化，得到在压缩感知框架下雷达成像聚焦的数学表达式；

(5)根据发射-接收-采样的三维形式将接收信号写成张量形式；

(6)对张量接收信号作高阶奇异值分解，得到多维线性测量结果；

(7)采用张量混合匹配追踪方法对步骤(6)得到的稀疏信号恢复；

(8)将恢复的向量按照预先划分好的网格进行矩阵化处理，得到最终MIMO雷达稀疏成像的结果；

(9)在色噪声情况下，划分两个子发射阵列，构造互协方差张量，通过高阶奇异值分解去除色噪声带来的不利影响。

所述步骤(5)的张量形式建立过程如下：

(5.1)获得单基地共址MIMO雷达空间谱域回波：

(5.2)划分成像网格点，得到离散稀疏信号模型；

且有

(5.3)将上述接收信号按照发射-接收-采样的三维信息写成张量形式

所述步骤(7)中所述采用张量混合匹配追踪方法对得到的稀疏信号恢复的步骤如下：

(7.1)初始化；首先定义支撑集

Λ_old＝max_ind(|σ_kron-omp|，K)，

其中σ_kron-omp＝kron-omp(Z，B₁，B₂，B₃，K)定义为标准kron-OMP方法的计算结果；残差初始化为

(7.2)支撑集扩充至2K个；

其中，

(7.3)支撑集更新；新的支撑集为

Λ_new＝max_ind(Z×₃B₁(Λ_temp)^T×₂B₂(Λ_temp)^T×₁B₃(Λ_temp)^T,K)；

(7.4)残差更新；

(7.5)迭代终止判断；通过迭代来持续更新残差和支撑集，当残差范数满足误差容限时，迭代停止，计算和输出σ；

所述的步骤(9)去除色噪声影响按照如下方法进行：

(9.1)划分发射阵列为两个子阵列，第一个子阵包含发射阵列的前M₁个天线，第二个子阵包含剩余的M₂＝M-M₁个天线；

(9.2)分别进行匹配滤波处理，得到

(9.3)将每个脉冲的匹配滤波输出堆栈成一个向量

(9.4)根据张量定义构建3阶张量

(9.5)根据(9.4)中的3阶张量，定义4阶协方差张量

(9.6)对协方差张量进行高阶奇异值分解去除色噪声影响

本发明的有益效果在于：

1.与DAS类成像方法相比，本发明克服了DAS类方法固有的低分辨率和高旁瓣的缺点。

2.与其他经典化压缩感知成像方法相比，本发明所提的THMP方法充分利用接收信号的张量特性进行稀疏信号恢复，避免了向量化操作带来的信号内在结构信息损失。

3.此外，本发明所提的THMP方法中的每一次索引选择的过程是利用OMP方法实现的，这种操作保证了在基信号选择时的正交性，也就能在字典矩阵具有傅里叶类似性质的时候可以区分相距很近的空间面元；与此同时，在THMP方法中存在的回溯选择操作与SP方法相同。这种操作的存在保证了THMP方法有能力剔除在前面的迭代过程中被选择的病态的索引，向支撑集中添加新的潜力高的索引。因此，THMP方法在理论上比OMP方法和SP方法的性能都要好。

附图说明

图1是采用HMP方法的MIMO雷达成像具体流程图；

图2是本发明的单基地MIMO雷达二维成像模型示意图；

图3是接收信号张量模型立体示意图；

图4是张量分解原理图；

图5是Kron-OMP方法的MIMO雷达的成像结果；

图6是NBOMP方法的MIMO雷达的成像结果；

图7是THMP方法的MIMO雷达成像结果；

图8是信号恢复的均方根误差与信噪比关系图；

图9是信号恢复的均方根误差与采样拍数关系图；

图10是色噪声情况下的信号恢复均方根误差与信噪比关系图；

图11是色噪声情况下稀疏信号恢复概率与信噪比关系图。

具体实施方式

下面结合附图对本发明作进一步描述。

本发明的目的在于克服上述技术的缺陷，提出一种基于张量混合匹配追踪的MIMO雷达稀疏成像方法。该方法利用MIMO雷达接收信号的多维结构，做高阶奇异值分解得到多维线性测量结果。然后在稀疏信号恢复框架下，将OMP方法和SP方法的优点结合起来，使得它在选择基信号时保证了正交性，而在支撑集更新时采用回溯策略。通过这种操作，所提方法在付出一定运算量的代价下可以保证很高的雷达图像重建分辨率，且不会出现伪影现象。

本发明成像方法主要包括以下几个方面：

1、推导出MIMO雷达的张量稀疏成像模型

如图2中所示，MIMO雷达由M个发射阵元和N个接收阵元组成，且发射阵元、接收阵元分布在二维平面上的同一基线上。MIMO雷达成像几何关系如图2所示。以成像场景中心为极坐标原点建立坐标系，则第m个发射阵元和第n个接收阵元可分别表示为和(R_Rx,m,φ_Rx,m)。其中是收发阵元与Y轴正向的夹角。设目标的第k个散射点的直角坐标为r_k＝(x_k，y_k)，散射系数是σ(r_k)。第m个发射天线到第k个散射点距离记为第n个接收天线到第k个散射点的距离为天线阵列基线到场景中心的距离是R₀。

第m个天线发射信号S_m(t)可表示为

S_m(t)＝p_m(t)exp(j2πf_ct) (1)

式中，p_m(t)是发射信号的归一化包络，f_c是载波频率。MIMO雷达发射的是相位编码正交信号，假设其具有理想的自相关特性和互相关特性。

设成像场景中共有K个散射点，则M个发射信号经过K个散射点反射到第n个接收阵元被接收的叠加回波为

其中τ_n，m(k)是第m个发射阵元到第k个散射点再到第n个接收阵元的整个辐射过程的路径延迟。根据远场假设，则有|r_k|＜＜R_Tx，m，|r_k|＜＜R_Rx，m，则第m个发射阵元到第k个散射点的距离和第k个散射点到第n个发射阵元的距离可近似为

其中I_Tx，m和I_Rx，n分别是第m个发射阵元和第n个接收阵元到成像场景中心的单位位置矢量，即

则延迟τ_n，m(k)可以近似为

其中，c是电磁波传播速度，上式中与R_Tx，m和R_Rx，n有关的项都属于固定已知项。

去载波后，通过相关处理器组(匹配滤波)，利用发射信号的正交性实现通道分离输出的第(n，m)个通道信号是

对其做傅里叶变换，并且代入路径延迟公式，得到频域形式的输出为

上式中，令

其中是MIMO雷达第(n，m)个观察通道的波数。

由于MIMO雷达存在许多不同收发组合，随着和的变化，它们将填充空间谱域一定范围的支撑区分布。因此，我们可以得到相位分集MIMO雷达在空间谱域的回波表达式

上式表明，经过匹配滤波后目标的散射系数与第(n，m)个通道在空间谱的回波满足傅里叶变换关系。设空间谱中每个通道均有q个采样样本，则式(9)可表示为向量形式

z_n，m＝[z_n，m(K_n，m(f₁))…z_n，m(K_n，m(f_q))]＝A_n，mσ (10)

其中z_n，m(K_n，m(f_i))是第(n，m)个通道的在空间谱域的第i个观测样本，而A_n，m是第(n，m)个通道的观测矩阵，σ∈C^K×1是K个散射点构成的向量，且有

公式(9)表明目标的后向散射系数(Radar Cross Section，RCS)和MIMO雷达接收信号在空间谱域的表达式是一对傅里叶变换对。假定在每一个信道有Q个采样点。在匹配滤波和傅里叶变换后，在每个采样时刻，数据都呈矩阵形式。传统的矩阵分析方法都会堆栈公式(9)中每一个采样点的数据来形成一个大的接收数据矩阵。这类方法忽视了发射阵列、接收阵了h采样序列的多维结构。根据张量的定义，多个采样点的接收数据可以写成张量形式。根据张量的矩阵展开定义，三阶张量数据可以展开为

由公式(12)可知，MIMO雷达张量接收数据的模-3展开的转置与接收信号的矩阵形式相同。具体的接收信号张量结构见图3。

2、在MIMO雷达张量回波模型的基础上，利用高阶奇异值分解得到MIMO雷达多维线性测量结果。

一个原始的张量信号Y能作高阶奇异值分解：

Y＝X×₁D₁×₂D₂×₃D₃ (13)

其中，X为分解核张量，D₁，D₂和D₃是分解的因子矩阵，具体张量分解结构如图4所示。

当因子矩阵满足正交性时，核张量可以由下式得到：

因此，我们可以选择合适的因子矩阵使得核张量越稀疏越好。

高阶奇异值分解表达式(14)可以通过堆栈张量写成向量形式：

其中x＝vec(X)，y＝vec(Y)。

对稀疏张量的测量过程可以通过多线性测量解决。在各个维度进行线性测量：

Z＝Y×₁Φ₁×₂Φ₂×₃Φ₃ (16)

Z＝X×₁Φ₁D₁×₂Φ₂D₂×₃Φ₃D₃ (17)

上式也能写成向量化矩阵的形式

其中，z＝vec(Z)，B_n＝Φ_nD_n(n＝1,2,3)。K是信号的稀疏度。

3、在MIMO雷达有噪回波模型的基础上，利用混合匹配追踪方法进行MIMO雷达成像重构。

采用标准Kron-OMP方法得到稀疏解的初值且确定初始支撑集

σ_kron-omp＝kron-omp(Z，B₁，B₂，B₃，K) (15)

上式是指在三个线性测量矩阵分别是B₁、B₂和B₃，稀疏度是K，测量向量是Z的情况下标准kron-OMP方法的输出。

初始支撑集为

Λ_old＝max_ind(|σ_kron-omp|，K) (16)

其中，max_ind(p,k)函数是指返回p中幅值最大的k个元素所对应的索引。

则残差初始化为

采用Kron-OMP方法对上式所求的残差进行处理，可以得到

利用上式将支撑集扩充至2K个

将原始测量信号向这2K个支撑集构成的子空间投影，可以得到更新的支撑集

利用上式更新残差

迭代更新残差值R和支撑集Λ，从而提高稀疏解的恢复精度。最终可以得到稀疏解

将支撑集扩充至2K个，新添加的支撑集索引是由将上一步Kron-OMP方法输出的结果最大元素的标号所组成的；然后将回波采样向量Z向子空间投影，计算投影系数，利用投影系数中前K项最大的值组成新的支撑集Λ_new，这就是支撑集更新的过程。这一步，与SP方法中支撑集更新类似，具有回溯特性。所提MIMO雷达稀疏信号恢复方法可以总结如方法1所示。

从上面方法的描述来看，THMP方法中的每一次索引选择的过程是利用OMP方法实现的，这种操作保证了在基信号选择时的正交性，也就能在字典矩阵具有傅里叶类似性质的时候可以区分相距很近的空间面元。与此同时，在HMP方法中存在的回溯选择操作与SP方法相同。这种操作的存在保证了THMP方法有能力剔除在前面的迭代过程中被选择的病态的索引，向支撑集中添加新的潜力高的索引。通过以上分析不难看出，THMP方法在理论上比OMP方法和SP方法的性能都要好。

4、在色噪声条件下，提出构造互协方差张量的方法来消除色噪声对所提方法产生的不良影响。

为了能够有效地利用匹配滤波器的正交性抑制MIMO雷达接收信号中的空域色噪声，首先将M个发射天线分为两个子阵，第一个子阵包含发射阵列的前M₁个天线，第二个子阵包含剩余的M₂＝M-M₁个天线。然后分别用前M₁个发射波形和后M₂个发射波形对接收信号进行匹配，则有

式中将每个脉冲的匹配滤波输出堆栈成一个向量，则有

式中D₁＝A₁⊙B，D₂＝A₂⊙B，本发明将考虑接收信号的固有多维结构特性，提出色噪声环境下基于张量互协方差张量分解的解决方法。

根据张量的概念可知，式(25)和式(26)中的接收数据分别可以构建3阶张量和且满足

根据式(27)中的两个3阶张量，定义一个4阶协方差张量如下所示

式中n,i＝1,...,N.q＝1,...,M₂.j＝1,...,M₁。根据空域色噪声的特性则空间色噪声矩阵和满足因此在式(28)中，由于采用了不同匹配滤波器输出噪声的正交特性，互协方差张量中空域色噪声的影响得到消除，即通过构造互协方差张量达到了消除空域噪声的目的。对互协方差张量进行高阶奇异值分解，则有

式中表示核张量，且核张量满足正交特性。均为酉矩阵。此后，就可用标准THMP方法进行数据恢复。

本发明的效果可通过以下仿真说明：

仿真条件与内容：

1、MIMO雷达对单点目标和多点目标成像性能

真实验中MIMO雷达的收发阵列都为均匀线阵，4个发射阵元在X轴上，坐标设为(0,4,8,12)×λ/2，接收阵元亦为4个，坐标为(0,1,2,3)×λ/2。发射波形采用循环方法(cyclic algorithm-new，CAN)设计的正交波形，每个发射波形包含码元数为100，载频为10GHz，带宽设为50MHz，相应的码元时宽是0。02μs，脉冲重复周期为6μs，采样周期等于码元时宽。

将成像区域网格化。设观测区域由50个距离单元组成，方位角范围为-80°～80°，角度单元设为5°。假设目标坐标分别为(40,0°)，(40,-20°)，(25,20°)，(5,-10°)，(5,10°)，(10,10°)，(15,-10°)，(15,10°)和(45,10°)。所有点目标的后向散射系数都设为1，设噪声为加性高斯白噪声，信噪比为30dB，图5、图6、图7分别给出了标准Kron-OMP方法，NBOMP方法以及本文所提的THMP方法的成像结果。

图5、图6和图7可以看出分别给出了在多个点目标存在的情况下不同方法的成像结果。可以看出，三种张量恢复方法都能反演出感兴趣场景的散射点分布情况，但是性能有所差异。为了更直观的比较各种方法的成像性能，图8给出了多个点目标时信号恢复的均方根误差与信噪比关系图。

从图5、图6、图7和图8可以看出，Kron-OMP方法和NBOMP在假设点目标位置都形成了峰值，表明Kron-OMP方法和NBOMP在多点聚焦成像的可用性。但是存在明显的伪影点，不利于目标的判读。这是由OMP方法对支撑集进行扩充而从不加以删除的策略造成的，而上述两种方法是OMP方法在张量情况下的直接推广。THMP方法在所有目标位置处形成峰值，且旁瓣水平和分辨率都较上述两种方法好。这与前面的分析是一致的。THMP方法中每一次支撑集选择过程是利用OMP方法实现的，这种操作保证了基信号选择的正交性，于此同时，在THMP方法迭代过程中与SP方法类似，可以去除在前面迭代过程中被选择的病态索引，因此HMP方法分辨率更高。

此外，图9给出了三种方法信号恢复的均方根误差与采样数关系图。可以看出，三种张量恢复方法的恢复误差都随着采样数增多而减少，这是由压缩感知理论所决定的。采样拍数越多，训练样本越多，得到的稀疏解精度越高。而本发明所提的THMP方法较另外两种方法性能更好，在同样精度下所需采样数更少，这意味着所提方法的优越性。

图10和图11体现了在色噪声情况下，本发明提出的处理方法的优越性。与之对比的是Chen所提出的截断奇异值分解方法(A new method for joint DOD and DOA estimation in bistatic MIMO radar.Signal Processing,2012,90:714-718)。其中，图10是色噪声情况下的信号恢复均方根误差与信噪比关系图，图11是色噪声情况下稀疏信号恢复概率与信噪比关系图。可以看出，相比于，本发明所提的张量高阶奇异值分解方法有更好的性能。这是因为MIMO雷达信号的张量结构可以提高对信号子空间分解的精确度以及本发明利用了空域色噪声的正交性来消除它的影响。