一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法

2021-07-05邓理康张双辉张刘永祥

雷达学报 2021年3期

邓理康张双辉张弛刘永祥

(国防科技大学电子科学学院长沙 410073)

1 引言

逆合成孔径雷达(Inverse Syntheic Aperture Radar,ISAR)可以生成目标的二维图像，通常用于目标识别和分类[1,2]。与二维图像相比，三维图像能够获得更多的目标信息，更便于目标识别。因此，近年来得到了越来越多学者的关注。文献[3]研究了干涉逆合成孔径雷达(Interferometric Inverse Synthetic Aperture Radar,InISAR)技术提取目标三维特征，但是需要面临复杂的运动补偿问题。文献[4]提出了多输入多输出(Multiple-Input Multiple-Output,MIMO)的三维成像方法，该方法采用十字交叉阵列分别获得2个维度的分辨，发射宽带雷达信号获得距离分辨，因此可以在单次快拍下获得三维图像从而避免了运动补偿。但是需要的阵元数量大，硬件成本高。文献[5]采用窄带MIMO雷达获得三维图像，阵型采用接收面阵和发射线阵的设计，发射窄带信号对目标进行单次快拍成像。但是为了获得高分辨的三维图像，需要的阵元数较多。文献[6]将ISAR技术与MIMO雷达相结合，在有限成像时间内，能够提高图像分辨率，并且用时间采样替代空间采样，降低了硬件成本。但是以上方法面临的共同特点就是需处理的数据量大，实现的硬件成本高，因此需要一种通过有限数据和硬件就能得到高分辨图像的方法。

压缩感知(Compressed Sensing,CS)理论[7–8]利用图像的稀疏性，用少量采样数据就能恢复出高分辨图像。因此，在图像处理领域获得了广泛应用。文献[9]利用克罗内克压缩感知(Kronecker Compressed Sensing,KCS)方法构建感知矩阵，并将多维信号转化为一维向量，从而实现了从多维CS到一维CS的转换。然而KCS存在的问题是测量矩阵和信号的维数过大，导致了存储和计算负担急剧增加，影响了KCS算法进一步应用。为了降低测量矩阵和信号维数，文献[10]首先将三维数据切片，再沿着一维叠加，将三维张量转换为二维矩阵。然后利用降维二维平滑l0范数(Dimension Reduction2D Smoothedl0-norm,DR-2D-SL0)重构方法求解稀疏优化问题，得到散射体的二维图像。最后，通过将得到的二维图像重新还原成三维张量来得到三维图像。文献[11,12]首先将信号表示为张量和矩阵的模态积形式，然后采用一种稀疏多维度平滑l0算法(MultiDimensional Smoothedl0-norm,MD-SL0)恢复三维图像，由于该算法能够直接处理三维图像，所以具有较低的计算复杂度和内存消耗，能够有效地进行三维重建。但是该方法需要调整的参数较多，并且在低信噪比条件下图像恢复效果一般。文献[13]将ADMM算法融入稀疏孔径成像ISAR自聚焦中，利用部分傅里叶矩阵字典将矩阵求逆转化为元素除法，并将二维ISAR图像的距离单元更新变为整体更新，在数秒内就能获得聚焦良好的图像。

受文献[13]的启发，本文基于ADMM算法提出了一种MIMO-ISAR成像方法。首先建立多维回波欠采样模型，然后利用ADMM算法解决欠采样数据的稀疏约束欠定问题。为了提高计算效率、降低存储消耗，该方法将高维的测量矩阵分解为张量的模态积，用张量元素除法替代计算效率低的矩阵求逆。仿真和实测数据都验证了所提方法的有效性。下面给出本文的符号定义，⊗表示两个矩阵的克罗内克积；×n(n=1,2,3)表示张量和矩阵的n模态积；张量用大写花体字母表示，如A。

2 信号模型

图1为成像场景图，其中X轴沿着收发阵列方向，选取一发射节点为坐标原点O，Y轴和Z轴分别按照右手定则确定。飞机以速度为v0作匀速直线运动，且速度方向与阵列方向相垂直。

图1 成像场景图Fig.1 Geometry of imaging

图1的线性收发阵元如图2所示，假设有M1个发射阵元和N1个接收阵元，其中发射阵元的间距为2N1d，接收阵元的间距为 2d。根据文献[14]，上述收发阵列可以等效成A=M1N1个等效收发阵元，且等效收发阵元成线性均匀排列，阵元间距为d。

图2 线性收发阵元的等效收发阵元示意图Fig.2 The equivalent transceiver array element for linear receiving and transmitting array elements

假设发射阵列发射步进频信号，发射频率fb=fc+(b −1)Δf，其中fc为发射信号中心频率，Δf为频率步长。经过信号分选后，第a个等效收发阵元在第p个快拍时刻的接收信号可以表示为

其中，a=1,2,···,A表示等效收发阵元序号，A为等效阵元总数；q=1,2,···,Q表示散射点序号，Q为散射点总数；p=1,2,···,P表示快拍序号，P为快拍数；c为光速；σq为第q个散射点的散射强度；表示第a个等效收发阵元在第p个脉冲时刻到第q个散射点的距离。

如图3所示，目标在ZOY面沿着Y轴方向以速度v0作匀速直线运动，假设初始快拍时刻，目标参考点位于O0，第p个快拍时刻参考点位于Op,φa ≈(a −1)d/R0+α0为第a个等效阵元与Z轴夹角，θp ≈ω(p −1)Tp为OOp与Z轴夹角。其中d为等效阵元间的距离，ω为目标运动等效转速，Tp为脉冲重复时间，R0为等效收发阵元中心到目标的距离，ω ≈v/R0,α0为第1个等效收发阵元与Z轴夹角。则在第p个快拍时刻，坐标系中q点坐标转变为坐标系中的坐标其中轴由O和O0的连线确定，轴和轴分别由右手定则确定；轴由第a个等效阵元和Op的连线确定，分别由右手定则确定。则经过平动补偿后，式(1)中的回波信号可以表示为

图3 目标移动示意图Fig.3 Motion of the target

其中，Rap表示第a个等效阵元到转动中心Op的距离。则

将式(4)代入式(2)可得

当发射信号带宽远小于中心频率fc时，可近似认为fb/c≈fc/c，假设α0≈0，在远场条件下且观测时间较短时，则d ≪R0,φa ≈0,θp ≈0。则式(5)可以近似为

则通过对x方向的相位求偏导数，可以得到沿x方向相位变化为

假设两个散射中心在x方向上的距离为 Δx，则总相位差可表示为

当满足ψx ≥2π时，两个散射点在x方向才能被分辨，所以x方向上的分辨率为

同理可得y方向的分辨率为

而z方向的分辨率为pz=c/(2Bw)，其中Bw为信号带宽。所以对a,b,p3个维度分别做傅里叶变换就能得到三维图像。所以信号的张量模型可以表示为

其中，×n(n=1,2,3)，表示张量和矩阵的n模态积，S ∈CA×P×B表示三维信号，X ∈CA×P×B表示三维图像。F1∈CA×A,F2∈CP×P,F3∈CB×B表示全傅里叶变换矩阵。

当阵元数目减少，或由于噪声或硬件原因导致完整回波出现缺失时(相当于对回波的3个维度做稀疏采样)，如果继续采用傅里叶变换将会导致图像质量严重下降。因此我们利用图像的稀疏性，引入压缩感知方法对回波信号进行处理，用较少的数据就能恢复图像。但是传统的压缩感知算法需要将三维数据展开成一维形式。将式(11)展开成一维形式为

其中，符号⊗表示克罗内克积，s=vec(S),x=vec(X)。假设三维信号的维度为60×60×60，当转化为一维信号后，感知矩阵的维度为21600×21600，这需要消耗巨大的内存和存储容量，限制了算法的进一步使用。所以针对上述问题，本文采用多维ADMM(MultiDimensional ADMM,MD-ADMM)算法对图像进行稀疏重构。

3 多维ADMM稀疏成像方法

3.1 多维ADMM算法

下面以ADMM算法为基础，推导MD-ADMM算法，首先推导算法的三维张量形式，再推广到多维情形。当完整回波出现缺失时，式(11)可以表示为

其中，S ∈CM×N×K表示降采样信号，X ∈CU×V×W为待恢复图像。设完整信号维数为U ×V ×W,M,N,K分别为对信号3个维度的采样数。F(1)∈CM×U,F(2)∈CN×V,F(3)∈CK×W，分别为3个维度的部分傅里叶变换矩阵。令F(1)=T1F1,F(2)=T2F2,F(3)=T3F3，其中F1,F2,F3分别为维数为U ×U,V ×V,W ×W的全傅里叶矩阵。T1∈CM×U,T2∈CN×V,T3∈CK×W为傅里叶采样矩阵。令G,H,J分别表示对张量S3个维度的采样序列，其中，G ∈[1,2,···,U]T,H ∈[1,2,···,V]T,J ∈[1,2,···,W]T其序列长度分别为M,N,K。则

将式(13)展开成一维形式为

其中，s=vec(S),x=vec(X)，其中vec(·)表示将张量向量化。在式(14)的信号模型中，直接从已知信号s中重构出x将有无穷多组解，属于病态问题，因此需要引入额外条件，压缩感知理论可以利用待恢复信号的稀疏性从无穷多组解中找出最稀疏的解。雷达图像通常由若干强散射点组成，因此待恢复信号x的稀疏性一般情况下是成立的。理想情况下通常用l0范数表示信号的稀疏性，但是基于l0范数最小化的稀疏恢复问题不易求解，属于NP难问题。因此，一般采用l1范数替代l0范数，并且l1范数最小化问题通常可以转化为凸问题。本文采用ADMM算法[15]解决基于l1范数的最小化问题，该方法可以将高维的测量矩阵分解为张量的模态积，用张量元素除法替代了计算效率低的矩阵求逆，提高了计算效率。因此式(14)中，基于l1范数最小化优化问题可以表示为

引入辅助变量z，则原基于l1范数的最小化问题可以等价于以下带等式约束的优化问题

其增广拉格朗日函数可以写为

其中，α为对偶变量，ρ为惩罚系数，ADMM算法将式(15)中的问题分解为如式(18)的3个子问题

其中，上标k表示第k次迭代，式(18)中前两个子问题可以通过对函数Lρ(x,z,α)中的x和z分别求导并令导数为0求得，经求解式(18)中x,z以及α的迭代式为

其中，ST为软阈值(Soft Threshold)函数，其表达式为ST(x,a)=(x/|x|)max(|x|−x,0)。将F(1)=T1·F1,F(2)=T2F2,F(3)=T3F3代入式(19)可得

将式(21)写成张量形式，可得

其中，1U×V×W表示元素全为1、维数为U ×V ×W的三维张量，⊘表示张量的元素除法，G表示对接收回波三维方向的采样，其值设为0或1，分别表示是否被采样到。同理式(19)中z(k+1)和α(k+1)的迭代式同样可以写成如式(23)和式(24)的张量形式

联合迭代式(22)—式(24)，就可得到图像X。初始参数设置如下：Z和A的初值设定为0，λ的值根据数据进行调整，ρ的值取1。

将三维形式进一步推广，假设多维张量的维数为U1×U2×···×Ui，则张量X(k+1)的更新表达式为

3.2 计算复杂度分析

MD-ADMM算法中图像X更新表达式(22)的计算复杂度为：O{[MNKU+UVNK+UVKW+2(U2VW+V2UW+W2UV)]假设经过L1次迭代后算法停止，则总的算法复杂度可以表示为：O{[MNKU+UVNK+UVKW+2(U2VW+V2UW+W2UV)]L1}。根据文献[10]，算法DR-2D-SL0的计算复杂度为O{(UVMNK+MNUVW+MNWK+UVKW)L2L3}，其中L2和L3分别为第1层循环和第2层循环的迭代次数。根据文献[5]，算法MD-SL0的计算复杂度为O[(MNKU+UVNK+UVKW+UVWM+MNVW+MNKW)L2L3]。在本实验中，迭代次数L2L3≫L1，且算法DR-2D-SL0的单次迭代计算复杂度最高，因此3种算法的计算复杂度排序为：MD-ADMM

4 实验仿真

4.1 仿真数据分析

实验仿真数据设置如下：目标飞行速度v0=200m/s，雷达距目标中心的距离R0=10000m，发射信号中心频率fc=10GHz，带宽Bw=150MHz，发射步进频信号个数B=60。设收发阵列为10发6收MIMO线阵，等效收发阵元个数A=60，等效收发阵元间距d=2.5m。脉冲重复频率PRF=80Hz，快拍数P=60。仿真中使用的点散射模型如图4所示。

图4 仿真目标三维散点图Fig.43 D scatter of simulation target

当回波数据完整时，对3个维度直接进行傅里叶变换后得到的图像如图5所示。由图5可知当回波数据完整时，直接对3个维度做傅里叶变换可以得到质量较高的三维图像。本文提取图5中的三维散点坐标，将散射点所在位置幅值设为1，图像中其余位置幅值设为0，形成参考三维图像H。下面对稀疏采样回波进行成像，为了获得三维稀疏回波，对回波进行稀疏采样，具体采样方式如图6所示。

图5 完整回波数据图像三视图Fig.5 Three views of image with the complete echo

图6 回波采样形式Fig.6 Undersampling masks of random sampling and block sampling

首先采用随机采样方式，对稀疏度(每个维度的稀疏度相同)分别为50.0%,33.3%,25.0%的回波进行三维成像处理，其中信号添加信噪比为20dB的高斯白噪声。分别采用RD,MD-SL0,DR-2D-SL0,MD-ADMM4种算法对目标进行三维成像(其中算法MD-SL0和算法DR-2D-SL0统称为SL0算法)，并采用三视图进行展示。在实验中本文调整参数让每一种算法都达到其最佳成像效果。图7—图9分别为稀疏度为50.0%,33.3%,25.0%的不同算法成像结果。由图7—图9可知，当回波稀疏时，由于受到旁瓣干扰，传统RD算法将会失效，得到的图像分辨率较低。由于利用了图像的稀疏性，采用了压缩感知方法，算法DR-2D-SL0,MD-SL0,MD-ADMM都得到了质量较高的图像。

图7 稀疏度为50.0%时图像Fig.7 Image when sparsity is50.0%

图8 稀疏度为33.3%时的图像Fig.8 Image when sparsity is33.3%

图9 稀疏度为25.0%时的图像Fig.9 Image when sparsity is25.0%

为了进一步定量比较4种算法，表1给出了随机稀疏采样条件下4种算法数值结果，其中包括图像熵、峰值信噪比(Peak Signal-to-Noise Ratio,PSNR)以及计算时间。在图像处理领域，图像熵和PSNR可以一定程度上反映图像质量。在三维图像中图像熵的定义为式(26)

其中，Sum(·)表示张量的所有元素之和，C=U,V,W分别为张量X的维数。假设对X进行了归一化，则X相对于参考图像H的均方误差可以表示为

其中，huvw为参考三维图像H的元素。定义PSNR为

由表1以及图7—图9可知，虽然RD算法计算效率最高，但是图像质量最差。由于参数设置一致，算法DR-2D-SL0和MD-SL0图像熵与PSNR相同，但是MD-SL0算法直接对三维张量进行处理，而DR-2D-SL0要将三维张量展开成二维矩阵，这增加了计算量，所以MD-SL0的计算时间更短，计算效率更高。与其余基于压缩感知算法相比，所提MDADMM算法在不同稀疏度下的图像熵最小，峰值信噪比最大，并且计算时间最短，这验证了所提算法的有效性。

表1 随机稀疏采样条件下数值结果Tab.1 Numerical results under random sparse sampling condition

接着采用块稀疏采样方式，对稀疏度(每个维度的稀疏度相同)为50.0%,33.3%,25.0%的回波进行成像处理。图10为采用不同算法对信噪比为20dB，稀疏度为50.0%的回波进行处理后得到的图像。表2为块稀疏采样条件下不同算法的数值结果，由图10和表2可知，所提MD-ADMM算法在不同稀疏度的块稀疏采样条件下得到的图像熵最小，峰值信噪比最大，并且计算时间最短，这验证了所提算法的有效性。

表2 块稀疏采样条件下数值结果Tab.2 Numerical results under block sparse sampling condition

图10 稀疏度为50.0%的块稀疏采样图像Fig.10 The image of a sparsity of50.0%by block sampling

下面比较不同算法在相同稀疏度(3个维度的稀疏度相同)，不同信噪比下的成像性能，其中回波信号采用随机采样方式，每个维度的稀疏度均为25.0%，并添加均值为0的高斯白噪声。图11—图13分别为在信噪比为–5,0,10dB条件下不同算法得到的目标三视图。表3为不同信噪比条件下的数值结果。由图11—图13和表3可知，在稀疏孔径条件下RD算法基本失效。在不同信噪比条件下，所提MD-ADMM算法图像熵最小，PSNR最大，并且计算时间最短，这证明了所提算法对噪声的鲁棒性最强。

表3 不同信噪比条件下数值结果Tab.3 Numerical results under different SNR conditions

图11 稀疏度为25.0%信噪比为–5dB的图像Fig.11 The image when sparsity is25.0%and SNR=–5dB

图12 稀疏度为25.0%信噪比为0dB图像Fig.12 The image when sparsity is25.0%and SNR=0dB

图13 稀疏度为25.0%信噪比为10dB的图像Fig.13 The image when sparsity is25.0%and SNR=10dB

4.2 实测数据分析

MIMO雷达实验系统目前还在搭建中，还存在发射信号带宽过窄、收发阵列同步性差等问题，导致MIMO-ISAR回波目前还难以获取，也是下一步要着重解决的问题。因此采用二维Yak-42飞机ISAR实测数据，以验证所提MD-ADMM算法在有限采样数据条件下的有效性。实验数据参数设置如下：雷达发射信号的载频为5520MHz，信号带宽为400MHz，快时间采样频率为10MHz，脉冲宽度为25.6µs，观测目标为Yak-42飞机。假设接收信号已经做了包络对齐、平动补偿以及自聚焦，共接收到256个脉冲，每个脉冲信号包含256个快时间采样。本文采样稀疏采样方式，抽取96个脉冲，以及128个快时间信号。图14为不同算法对二维稀疏信号成像结果。图15为图14信号中添加0dB的高斯白噪声的结果。表4为不同算法在不同信噪比条件下的数值结果。由图14、图15以及表4可知，RD算法基本失效，算法MD-SL0,MD-ADMM在二维稀疏采样条件下依然能够得到清晰图像，但相比MDSL0算法，MD-ADMM算法得到的图像熵最小，峰值信噪比最大，并且计算时间最短，进一步验证了所提算法的有效性。