高效的MIMO雷达运动目标三维成像方法

王伟，胡子英，岳佳男

（哈尔滨工程大学自动化学院，黑龙江 哈尔滨 150001）

1 引言

多输入多输出（MIMO, multiple-input multiple-output）雷达是近年来提出的一类新型雷达。区别于其他工作方式，MIMO雷达在发射端发射正交信号，在接收端通过匹配滤波实现信号的通道分离，这样的方式使它拥有更多的虚拟孔径和更高的阵列自由度，极大地提高了雷达成像的性能[1-2]。为了更好地呈现实际目标的三维结构及特征，近年来，三维成像技术开始在MIMO雷达中得到研究及发展，现阶段主要从优化阵元的排布方式[3]、提高成像精度和效率[4-5]等方面展开。但是这些研究主要围绕空间中的静态目标，在实际应用中，目标运动因多普勒频移带来的影响不可忽视。

首先，多普勒频移的产生会使回波信号的中心频率发生改变，导致估计结果与实际情况存在较大偏差。其次，对多普勒频率的采样使测量矩阵各列间的互相干性不再满足Welch界限制条件[6]，稀疏重构的恢复概率难以得到保证。文献[7]从速度带来的频移影响考虑，用迭代自适应估计方法（IAA,iterative adaptive approach）对距离-角度-多普勒的三维特征参数进行估计，但IAA的高运算复杂度使其难以应用于三维成像中。文献[8]在回波中考虑了速度带来的多普勒影响，并用变分贝叶斯的方法对对速度参数进行精确的估计，但多层贝叶斯模型使对速度求解和稀疏成像的迭代过程更加复杂。

对于多普勒频率采样带来的强互相干性问题，文献 [9-11]分析了各列之间的互相干性[12]对测量矩阵性能的影响，并以Welch界判定准则构造凸优化问题，找到最优投影矩阵，使测量矩阵的性能逼近最佳。文献[13]则在此基础上利用信号的先验信息简化了该凸优化问题的处理过程，预先给定加权矩阵，提高了算法效率。但是在这些方法中，投影矩阵及原始测量矩阵均参与运算。使整个优化过程尤为复杂，并且MIMO雷达稀疏成像中需要将区域每点的散射系数张成一个向量来处理，使整个模型的数据结构变得极大[14]，因此，运动目标三维成像亟需解决运算效率过低的问题。针对提升算法运算效率，文献[15]利用矩阵完备性（MC, matrix completion）理论构造具有低秩和均匀随机采样性质的矩阵，来代替测量矩阵实现三维成像，成像效果良好，但对运算效率提升的效果一般。文献[16]中将所得回波模型用基于克罗内克压缩感知（KCS,Kronecker compressive sensing）的方式表示，并在算法的初始化阶段利用匹配滤波得到的结果提取目标大致范围，借此来优化测量矩阵，以降低 KCS成像的运算复杂度。文献[17]则在初始化目标区域的基础上实现了三维成像的张量分解，利用此结论为Tucker分解提供了先验，在保证其唯一性条件的前提下提升了运算效率，但文献[16-17]中的方法普遍受噪声的影响比较大，不适合在偏低的信噪比下工作。

本文针对 MIMO雷达远场条件下的运动目标提出一种高效的三维成像方法。以多普勒频率作为一个新维度，在进行对目标区域的稀疏反演前，依据各维度的回波依次完成对目标的搜索，得到空间目标在各个维度的分布信息，并根据索引重构一个新的低维测量矩阵，缩小目标区域范围，有效地提高了运算效率。同时分析多普勒频率采样对测量矩阵相干性的影响，并预先构造一个投影加权矩阵，应用高精度贝叶斯方法实现对该投影矩阵的优化，得到目标沿多普勒维度的精确分布信息。仿真结果表明，所提方法可以极大地提升运算效率，实现高效成像，具有精确稳定的成像性能。

2 MIMO雷达运动目标回波模型

MIMO雷达发射阵列是由M个阵元构成的均匀线性阵列，以Ti(i=1,…,M)表示第i个阵元，且在该发射阵元采用的是Hadamard正交相位编码信号s(t) =p(t) ej2πfct+ϕi，其中，p(t)为信号的复包络，iiifc为发射信号载波频率，ϕi是信号的相位。此类编码信号通过调整发射信号的相位ϕi来保证各发射信号间的正交性。接收阵列为由NR×NC个阵元构成的均匀分布二维阵列[16]，以Rj′l(j′ = 1 ,… ,NR；l=1,…,NC)表示接收阵元位于接收阵列的第j′行第l列。假设空间中的目标服从Sweiling II起伏模型，即对空间中K个远场运动目标点，在一个脉冲周期内散射系数不变，并以目标对应的多普勒频率fdk作为新的变量来表示目标运动带来的影响[18-19]，则第i个发射信号经过空间中K个运动目标点的反射后，由Rj′l接收到的信号为

其中，σk为目标点k的散射系数，τk,i,j’,l为信号经由第k个散射点的时延，由目标及阵元位置决定。经过去载波处理后可得接收端回波信号，如式(2)所示。

其中，q=1,…,Q表示第q次采样，tq为对应采样时间。按照压缩感知成像理论，对目标区域划分网格点，沿目标的3个维度划分的网格点数分别为U1、以表示网格点处的多普勒频率，则式(2)的回波实际应为

根据式(3)的模型进行稀疏成像，需依据式中的旋转矢量构造测量矩阵，就大多数研究而言，应用于CS成像方法中的测量矩阵是可以预先设计的，但式(3)中多普勒频率的存在使测量矩阵中出现了未知项，无法实现对测量矩阵的精确构建，且矩阵中的不确定性因素会对稀疏目标的恢复产生严重的影响，因此需要在对目标的三维坐标进行估计的同时对目标的多普勒频率进行精确的估计，故本文将目标的多普勒频率作为X、Y、Z坐标轴外的另一变量进行处理，且对式(2)进一步处理，得

其中，Tik表示Ti发射阵元到第k个目标点的距离，Rj'lk表示第k个目标点到Rj'l接收阵元的距离。由于实际中目标速度远小于信号传播速度，即所以可近似处理因此有

在成像区域中选定点O作为参考中心，且令TiO和Rj’lO分别表示O点到发射阵列、接收阵元的距离，则由O点反射的回波信号可作为参考信号实现对目标回波信号的补偿，即

如图1所示，远场目标区域中共有K个目标点，以O点为参考中心，在MIMO雷达收发阵列中，以T0和R00分别作为发射阵列、接收阵列的参考阵元。

图1 MIMO雷达远场成像系统几何模型

在文献[16,20]中对式(6)中所示的距离项做了等效近似处理，如式(7)所示。

综上，可以构造运动目标稀疏回波模型为

其中，有

且对测量矩阵有

3 MIMO雷达运动目标的高效三维成像

传统的 CS稀疏成像，需在各维度预先划分网格点，并在此基础上构成测量矩阵。通过稀疏恢复得到的稀疏解中包含了目标在区域内的分布信息，求解散射稀疏向量的稀疏恢复过程依赖测量矩阵每列与每个网格点的对应关系。对于稀疏成像而言，空间中大部分的位置散射强度认定为0，说明在稀疏恢复的过程中，测量矩阵的大部分列并不对应目标分布，却仍然要进行算法运算，这增加了算法的运算复杂度。在式(12)～式(14)模型中，上述无关列的无用运算量问题更明显，因此本文选择对测量矩阵进行处理，从中仅选取稀疏相关列，将大部分的无关列剔除，那么稀疏求解的运算效率就会大大提升。同时，由于测量矩阵的互相干特性主要由多普勒维度的采样所影响，基于本文方法可以单独在多普勒维度实现对测量矩阵的优化，提升整体算法的优化效率。

3.1 目标区域优化搜索

首先要完成的是得到所需要的相关列的索引，这些索引对应空间中目标点的位置，因此可以利用回波从各个维度搜寻目标信息，得到在每个维度目标的粗略分布情况，以构建索引支撑集。那么各维度都可以得到一个典型的CS模型，如式(15)所示。

如此构成了各维度的稀疏恢复问题，以在X维的搜索为例，上述对应X维的模型可以转化为对如式(18)所示的凸优化问题的求解。

各种方法（如正交匹配追踪（OMP, orthogonal matching pursuit）算法、稀疏贝叶斯学习（SBL,sparse bayesian learning）方法等）均可以实现对以上凸优化问题的求解，得到目标区域沿X维的分布信息，在划分的X维网格x1,x2,…,xUX中，可以得到对应目标在X维的位置索引构成的集合{ΛX|Λx1,Λx2,…,ΛxPx,Px≤K}，Px≤K是考虑到了空间中不同目标点可能会有相同的X维坐标，较易引起搜索过程中目标点的丢失，这Px个位置包含了空间中所有的K个目标点的X坐标，如此便完成利用接收回波在X维度上对目标的搜寻工作。

同理，经过类似的处理可以得到其余3个维度上的坐标索引集合分别为{ΛY|Λy1,Λy2,…,ΛyPy，Py≤K}、{ΛZ|Λz1,Λz2,…,ΛzPz，Pz≤K}和{Λf|Λf1,Λf2,…,ΛfPf，Pf≤K}，将所得的 4 个集合进行外积处理，可以得到Px×Py×Pz×Pf个组合指向成像区域内的Px×Py×Pz×Pf个点，这些目标点的集合Λ中必然包含了所有的真实目标点，即Px×Py×Pz×Pf≥K，于是可将此集合作为新支撑集的索引构造测量矩阵ψ∈ CMNRNC Q×(Px×Py×Pz×Pf)， 式(12)中的压缩感知模型可转换为如式(19)所示的简化模型。

而P≪U，故采用常用的稀疏恢复算法时能实现良好的提升运算效率的效果。此外，由于各维度的搜索过程是依靠稀疏恢复算法完成的，而大多数算法可以实现良好的抗噪性，因此本文所提方法在低信噪比下也能正常工作。但需要注意两点，一是各维的CS模型需要保证存在唯一解；二是由式(8)可知，经过上述过程得到的稀疏解中，多普勒频率是目标在成像区域真实的多普勒频率，但其余坐标并不是目标在成像区域的真正坐标值，将其余坐标转换为实际坐标还需进行相应处理。基于式(8)和式(9)，令则目标其余坐标的真实坐标和得到的坐标之间的关系如式(21)所示。

因此，根据式(21)可得目标其余坐标在成像区域的目标坐标值。

经过上述处理后，得到了目标的三维坐标及多普勒频率，即实现了对运动目标三维成像算法运算效率的提升，但多普勒频率的引入会对测量矩阵的性能造成影响，对此需要分析测量矩阵性能的变化，并选取合适的方法来消除此影响，改善测量矩阵的性能，获得良好的稀疏恢复效果。

3.2 互相干性优化

本节对测量矩阵的性能进行优化，以实现良好的稀疏恢复。在常用的压缩感知模型中，常用约束等距特性（RIP, restricted isometry property）[12]来衡量性能，但是其在实际中的应用比较困难，因此可以用运算较便捷的互相干性来代替 RIP作为判断的准则。对一个M′×N′维的测量矩阵Θ，可以用Gram矩阵G=ΘHΘ来直观描述各列之间的互相干性，在对测量矩阵Θ的每列进行归一化处理后，最理想的情况下G应为单位矩阵，即只有对角线元素为1，其他位置元素均为0，代表各列向量之间是正交的。实际中这种完全理想的矩阵是极难得到的，因此能实现稀疏重构的较理想矩阵G[21]应该满足对角线元素为1，并且对角线外其他位置的元素必须不超过Welch界[8]的限制，然而在式(19)的压缩感知模型难以保证其性能。在X、Y、Z三维中，对目标信息的搜索没有问题，而影响整体测量矩阵性能的主要因素出现在多普勒维度的搜索中，由式(16)可知，多普勒维度上Q个采样点对应的采样时间而远场条件下采样时间的选取存在一个弊端，如图 2所示。

图2 远场采样等效距离

图2中，O为成像区域的中心，l1和l2表示成像区域与雷达间的最大距离和最小距离，R0表示成像中心到雷达的参考距离，Δl1、Δl2分别表示l1、l2与R0的距离差。信号发出后，在远场条件下接收端最早到最晚接收到信号的时间差为

由此可见，采样时间的取值与目标区域范围大小及信号传输速度有关，但要保证探测的精度，目标范围就不宜较大且实际中c认定为光速，这就导致Δt极小，而基于Q次采样均是在该时间内完成的，所以各采样时间取值均过小且在同一量级，使测量矩阵Af中的采样数据比较相近，致其相关性较差，难以保证稀疏重构的可靠性。

其中，g表示格拉姆矩阵。

首先，Af先经过列归一化处理，然后求解式(23)得到的Φ、g的精确联合估计，使测量矩阵的性能最佳。在解决凸优化问题的常用算法中，Bayes方法依据先验分布信息及后验概率密度实现对参数的精确估计，可获得较高的稀疏重构精度并且对噪声的抑制效果较好。在本文中，依据Gram矩阵的理想单位阵形式及投影矩阵的结构，确定复高斯分布先验信息，并且在多普勒维度的模型中，测量矩阵的维度较小，在处理过程中对应的乘法或取逆运算的运算量较小。因此，与文献[9-11]中在 Welch界准则下直接求解投影矩阵的方法相比，应用Bayes方法可以实现精度和效率二者之间的均衡，依靠复高斯分布先验信息可以获得更高的精度，同时本文对测量矩阵的降维处理也降低了运算的复杂度，所以本文从构造Bayes模型的角度考虑优化测量矩阵。在第s次迭代中，有

在上述目标函数下，需要做的就是在每次迭代中 Gram矩阵已确定的情况下获取对投影矩阵的最优解的精确求取。E中元素服从方差为γ的零均值复高斯分布，即Ea,b～CN(0,γI)，对应的概率密度函数为

投影矩阵Φ中元素也服从方差为β的零均值复高斯分布，即Φμ,κ～ CN(0,βI)，对应概率密度函数为

由此可得各参数的联合概率密度函数为

其中，复高斯分布的方差在[0,∞)范围内，即P(γ)∝1且P(β)∝1，因此式(26)可继续化简为

对式(28)取负对数处理可以得到各参数的目标函数，如式(29)所示。

经过推导，式(23)所示的原始问题转换为目标函数L(Φ,γ,β)的极小值求解问题，因此需将L(Φ,γ,β)对各参数求偏导，且其联合估计是以交替更新的方式实现的，以s表示第s次迭代，可得

3.3 运算复杂度分析

针对式(12)中的模型，若直接采用OMP算法，复杂度为O(I× (MNRNCQ) ×(UXUYUZUF))；若采用贝叶斯压缩感知（BCS, Bayesian compressive sensing）方法，复杂度为O(I2 ×(UXUYUZ UF))，其中I是算法运行的迭代次数。采用文献[16]的 DR-KCS方法时，其运算复杂度的高低与匹配滤波初始化得到的目标范围有关，DR-KCS应用于OMP算法时的复杂度为O(I× (MNRNC Q) ×P0)，应用于贝叶斯压缩感知方法时复杂度为O(I2×P0) ，其中P0与P的物理意义一样，都是经优化缩减后的目标范围，在高/低信噪比下P0的大小会有变化，进而会对运算复杂度产生影响。采用降维 Tucker分解[17]时，其复杂度近似为O(I× (MNRNCQ) ×(UXUYUZUF))。

本文方法的运算量分为两部分，一部分是在各维的搜索过程中，以X维为例，若采用OMP算法，复杂度为O(I×M×UX)，若采用贝叶斯压缩感知方法，其复杂度为O(I2×UX)，由此可见这个过程中的运算量极小，远低于下一过程，可忽略不计；另一部分是在式(19)的稀疏重构过程中，若采用OMP算法，复杂度为O(I× (MNRNCQ) ×P)，若采用贝叶斯压缩感知方法，复杂度为O(I2×P)。

本文方法与DR-KCS方法的复杂度相近，而与其他算法相比，由于P≪U，本文方法的运算量比较小，对运算效率的提升比较明显。

4 仿真实验与分析

本节的仿真旨在证明本文方法可以有效地降低高运算量对成像性能的影响，实现高效的成像。实验参数设置如下：发射阵元个数M=15，接收阵元行数与列数分别为NR=15，NC=15，发射阵元间距dt=2 m，接收阵元间距dr=1 m，信号载波频率fc=10 GHz，信号带宽B=50 MHz，发射脉冲宽度Tr=2 μs，采样次数Q=10，成像区域参考中心到收发阵列的距离为 10 km，并且依据文献[21]中对成像区域范围的限定，取X、Y、Z维的最大成像范围分别为-50～50 m，-50～50 m和 100～200 m，在该区域内设定由5个散射强度为1.0的强散射点构成，对应点的X、Y、Z三维坐标分别为(-20,10,150)、(0,-20,140)、(-30,0,130)、(20,-10,120)和(10,40,180)，5个目标点对应的多普勒频率分别为 0、666.7 Hz、1 166.7 Hz、1 500 Hz和 166.7 Hz。

图3描述的是在对测量矩阵优化前后各列之间的互相干性，其中，横坐标是各列互相干性的绝对值，互相干性绝对值为0表示两列完全不相干，其中互相干性绝对值最大值为 1，纵坐标表示与某一互相干性的值对应的列数。从图 3(a)中可以看出，在测量矩阵优化前，互相干性的绝对值为0～1的分布比较杂乱，并且相对来说，绝对值比较大，说明此时测量矩阵的互相干特性比较差，这是由于未对多普勒维度进行针对性优化导致的。而在图3(b)中，大多数列之间的相干性的绝对值为0，即大多数列之间互不相干，说明经过所提贝叶斯方法的优化处理后，测量矩阵的性能有了较大的提升，更有利于保证良好的稀疏恢复。

图4和图5是当信噪比SNR=30 dB时，本文所提方法的成像仿真结果，四维信息难以从图像中直接呈现，因此直接选取如图 4所示的XYFd、XZFd三维图像及如图5所示的XY、XZ、YZ这3个平面的投影来观察成像效果。本文所提方法可以很精确地估计散射点的位置和幅值，有很高的恢复概率。

图3 优化前后测量矩阵互相干性对比

图4 三维成像结果

图5 二维投影成像结果

图6展示了SNR=30 dB且目标个数为5时本文方法的估计结果与迭代次数的关系，其中实验参数与上述几次实验相同，针对性地从成像RMSE（root mean squared error）随迭代次数的变化关系来分析本文方法的收敛性，并且从图6中可以看出，本文方法在迭代5次左右基本完成收敛且误差逼近于0，说明有较好的收敛特性。

图6 成像误差随迭代次数变化关系

图 7(a)是描述成像误差与 SNR之间关系的曲线，其中本文方法应用BCS方法进行稀疏恢复。在SNR从-20～20 dB的范围内，对本文方法与OMP、MAP这 2种传统成像算法的稀疏成像性能做了比较，可以看出本文方法的成像性能稍有下降，这是因为重构测量矩阵时舍弃了大量的数据量，成像精度势必会受到影响。图7(b)是在SNR=30 dB下比较本文方法与MAP方法运行效率与稀疏度变化的关系，明显可看出本文方法效率极高，会节省很多的算法运行时间。

图7(c)是SNR在-20～20 dB的范围内对本文方法与 Matrix Completion (MC)[15]及降维压缩感知（DR-KCS, dimension-reduction kronecker compressive sensing）[16]这2种运算效率提升算法的成像性能对比，其中，本文方法与 DR-KCS方法均应用BCS方法进行稀疏恢复。从图7(c)中可看出，本文方法的精度更高，这是由于在各维的搜索本身就是稀疏恢复的过程，对噪声有一定的抑制作用，最后再经对式(19)的稀疏重构，因此会有更准确的稀疏成像效果。

图7(d)和图7(e)分别在不同的稀疏度和信噪比下比较了本文方法与DR-KCS方法的运行效率，图7(d)中SNR=30 dB，图7 (e)中目标数为5个。从图7(d)和图7(e)中可以看出，本文方法在目标个数较少的应用背景下有良好的处理效率，并且基于其各维搜索时的噪声抑制过程，使其在低信噪比下的性能明显优于DR-KCS方法。但目标个数多时算法效率会有所下降，这是由于新构造的测量矩阵的维度与目标个数成倍数关系，在多目标情况下测量矩阵列数会增大，因此使效率下降。

综合来看，由图7(a)和图7(b)对比可知，本文方法对一般的成像方法而言，运算效率提升的效果极为明显，虽精度稍有下降但可以满足实际的应用。在图7 (c)的对比中，本文方法与各效率提升算法相比又有更高的成像精度。在图7(d)和图7(e)的对比中，能够展现本文方法的运算高效性，尤其在目标个数较少及低信噪比下可以较好地实现对算法运算效率提升的目的。因此从整体来看，本文方法可以达到本文的设计初衷，即在尽量保证成像效果的条件下尽可能地提升算法的运算效率。

5 结束语

图7 本文方法与各类算法的性能对比结果

本文主要的贡献在于基于 MIMO雷达提出了一种针对远场运动目标的高效三维成像方法。从对测量矩阵处理的角度出发，依据回波模型依次在各维完成对目标的搜索，以初步缩小目标区域范围，实现对算法运算效率的提升。同时分析多普勒采样对测量矩阵相干性的影响，构造合适的投影矩阵。实现对整体测量矩阵互相干特性的优化。仿真结果表明，该方法可以在保证成像效果的条件下提升算法的运算效率，具有精确稳定的成像性能。