一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法邓理康      张双辉*      张  弛      刘永祥(国防科技大学电子科学学院   长沙   410073)摘   要：基于傅里叶变换的传统逆合成孔径雷达(ISAR)成像方法存在数据存储量大、数据采集时间长的问题。压缩感知(CS)理论利用图像的稀疏性，可以利用有限的数据恢复图像，这极大降低了数据采集成本。但对于多维数据，传统压缩感知方法要将多维数据转化成一维向量，这造成了很大存储和计算负担。因此，该文提出一种基于多维度-交替方向乘子法(MD-ADMM)的多输入多输出-逆合成孔径雷达(MIMO-ISAR)成像快速稀疏重建方法。首先建立基于张量信号的压缩感知模型，然后用ADMM算法对模型进行优化，将测量矩阵分解为张量模态积，用张量元素除法替代矩阵求逆，显著减少所需的内存和计算负担。该方法只需少量的数据采样，就能实现快速成像。与其他基于张量的压缩感知方法相比，该方法具有鲁棒性强、图像质量好、计算效率高的优点。仿真和实测数据验证了该方法的有效性。关键词：多维度-交替方向乘子法；压缩感知；多输入多输出-逆合成孔径雷达中图分类号：TN957.51 文献标识码：A 文章编号：2095-283X(2021)x-0001-16DOI: 10.12000/JR20132引用格式：邓理康, 张双辉, 张弛, 等. 一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法[J].雷达学报, 待出版. doi: 10.12000/JR20132.Reference format: DENG Likang, ZHANG Shuanghui, ZHANG Chi, et al. A multiple-input multiple-outputinverse synthetic aperture radar imaging method based on multidimensional alternating direction method ofmultipliers[J]. Journal of Radars, in press. doi: 10.12000/JR20132.A Multiple-Input Multiple-Output Inverse Synthetic Aperture RadarImaging Method Based on Multidimensional AlternatingDirection Method of MultipliersDENG Likang      ZHANG Shuanghui*      ZHANG Chi      LIU Yongxiang(College of Electronic Science and Technology, National University of Defense Technology,Changsha 410073, China)Abstract: The disadvantages of the traditional Inverse Synthetic Aperture Radar (ISAR) imaging method basedon Fourier transform include large data storage and long collection time. The Compressive Sensing (CS) theorycan use limited data to restore an image with the sparsity of the image, reducing the cost of data collection.However for multidimensional data, the traditional compressive sensing methods need to convert three-dimensional data into a one-dimensional vector, causing the storage and calculation burden. Therefore, thisstudy proposes a fast MultiDimensional Alternating Direction Method of Multipliers ((MD-ADMM)) sparsereconstruction method for Multiple-Input Multiple-Output ISAR (MIMO-ISAR) imaging. The CS model basedon the tensor signal was established, and the model with the ADMM algorithm was optimized. The measuredmatrix is decomposed into a tensor modal product, and matrix inversion is replaced by tensor element division,significantly reducing memory consumption and computational burden. Fast ISAR imaging can be achieved bya small amount of data sampling by the proposed method. Compared with other tensor compressed sensingmethods, this method has the advantages of stronger robustness, higher image quality, and computationalefficiency. The effectiveness of the proposed method can be invalidated by simulated and measured data.  收稿日期：2020-10-19；改回日期：2021-01-27*通信作者： 张双辉 shzhang3@126.com        *Corresponding Author: ZHANG Shuanghui, shzhang3@126.com基金项目：国家自然科学基金(61801484, 61921001)Foundation Items: The National Natural Science Foundation of China (61801484, 61921001)责任主编：张弓     Corresponding Editor: ZHANG Gong第10卷第x期 雷  达  学  报 Vol. 10No. x2021年1月 Journal of Radars Jan. 2021网络首发时间：2021-01-20 15:49:43网络首发地址：https://kns.cnki.net/kcms/detail/10.1030.TN.20210120.1419.002.htmlKey words: MultiDimensional Alternating Direction Method of Multipliers (MD-ADMM); Compressive Sensing(CS); Multiple-Input Multiple-Output Inverse Synthetic Aperture Radar (MIMO-ISAR)1    引言逆合成孔径雷达(Inverse Syntheic ApertureRadar, ISAR)可以生成目标的二维图像，通常用于目标识别和分类[1,2]。与二维图像相比，三维图像能够获得更多的目标信息，更便于目标识别。因此，近年来得到了越来越多学者的关注。文献[3]研究了干涉逆合成孔径雷达(Interferometric InverseSynthetic Aperture Radar, InISAR)技术提取目标三维特征，但是需要面临复杂的运动补偿问题。文献[4]提出了多输入多输出(Multiple-Input Mul-tiple-Output, MIMO)的三维成像方法，该方法采用十字交叉阵列分别获得2个维度的分辨，发射宽带雷达信号获得距离分辨，因此可以在单次快拍下获得三维图像从而避免了运动补偿。但是需要的阵元数量大，硬件成本高。文献[5]采用窄带MIMO雷达获得三维图像，阵型采用接收面阵和发射线阵的设计，发射窄带信号对目标进行单次快拍成像。但是为了获得高分辨的三维图像，需要的阵元数较多。文献[6]将ISAR技术与MIMO雷达相结合，在有限成像时间内，能够提高图像分辨率，并且用时间采样替代空间采样，降低了硬件成本。但是以上方法面临的共同特点就是需处理的数据量大，实现的硬件成本高，因此需要一种通过有限数据和硬件就能得到高分辨图像的方法。压缩感知(Compressed Sensing, CS)理论[7–8]利用图像的稀疏性，用少量采样数据就能恢复出高分辨图像。因此，在图像处理领域获得了广泛应用。文献[9]利用克罗内克压缩感知(KroneckerCompressed Sensing, KCS)方法构建感知矩阵，并将多维信号转化为一维向量，从而实现了从多维CS到一维CS的转换。然而KCS存在的问题是测量矩阵和信号的维数过大，导致了存储和计算负担急剧增加，影响了KCS算法进一步应用。为了降低测量矩阵和信号维数，文献[10]首先将三维数据切片，再沿着一维叠加，将三维张量转换为二维矩阵。然后利用降维二维平滑l0范数(Dimension Re-duction 2D Smoothed l0-norm, DR-2D-SL0)重构方法求解稀疏优化问题，得到散射体的二维图像。最后，通过将得到的二维图像重新还原成三维张量来得到三维图像。文献[11,12]首先将信号表示为张量和矩阵的模态积形式，然后采用一种稀疏多维度平滑l0算法(MultiDimensional Smoothed l0-norm,MD-SL0)恢复三维图像，由于该算法能够直接处理三维图像，所以具有较低的计算复杂度和内存消耗，能够有效地进行三维重建。但是该方法需要调整的参数较多，并且在低信噪比条件下图像恢复效果一般。文献[13]将ADMM算法融入稀疏孔径成像ISAR自聚焦中，利用部分傅里叶矩阵字典将矩阵求逆转化为元素除法，并将二维ISAR图像的距离单元更新变为整体更新，在数秒内就能获得聚焦良好的图像。n(n= 1;2;3)nA受文献[13]的启发，本文基于ADMM算法提出了一种MIMO-ISAR成像方法。首先建立多维回波欠采样模型，然后利用ADMM算法解决欠采样数据的稀疏约束欠定问题。为了提高计算效率、降低存储消耗，该方法将高维的测量矩阵分解为张量的模态积，用张量元素除法替代计算效率低的矩阵求逆。仿真和实测数据都验证了所提方法的有效性。下面给出本文的符号定义，表示两个矩阵的克罗内克积； 表示张量和矩阵的模态积；张量用大写花体字母表示，如。2    信号模型Ov0图1为成像场景图，其中X轴沿着收发阵列方向，选取一发射节点为坐标原点，Y 轴和Z 轴分别按照右手定则确定。飞机以速度为作匀速直线运动，且速度方向与阵列方向相垂直。 O发射节点接收节点(xTm,yTm)(xRn,yRn)ZYv0X 图 1 成像场景图Fig. 1  Geometry of imaging M1 N12N1d 2d A=M1N1 d图1的线性收发阵元如图2所示，假设有 个发射阵元和 个接收阵元，其中发射阵元的间距为 ，接收阵元的间距为 。根据文献[14]，上述收发阵列可以等效成 个等效收发阵元，且等效收发阵元成线性均匀排列，阵元间距为。fb = fc+ (b 1)∆f fc ∆f假设发射阵列发射步进频信号，发射频率，其中为发射信号中心频率， 为2 雷  达  学  报 第10卷ap频率步长。经过信号分选后，第个等效收发阵元在第个快拍时刻的接收信号可以表示为sap(fb) =Q ∑q=1qexp[ j4πfbcRqap](1)a= 1;2;


;A q= 1;2;


;Q Qp= 1;2;


;PPcq qRqap a其中， 表示等效收发阵元序号，A为等效阵元总数； 表示散射点序号，为散射点总数； 表示快拍序号，为快拍数；为光速； 为第个散射点的散射强度； 表示第个等效收发阵元在第p个脉冲时刻到第q个散射点的距离。v0 O0Op φa (a 1)d/R0 +

0 ap!(p 1)Tp OOp d!Tp R0!v/R0
0pO0-^x0^y0^z0q(^ xq0; ^yq0; ^zq0) Op-^xap^yap^zap(^ xqap; ^yqap; ^zqap)^z0 OO0^x0^y0^zap aOp^xap^yap如图3所示，目标在ZOY 面沿着Y 轴方向以速度作匀速直线运动，假设初始快拍时刻，目标参考点位于 ，第p个快拍时刻参考点位于 , 为第个等效阵元与Z 轴夹角，为 与Z 轴夹角。其中为等效阵元间的距离，为目标运动等效转速，为脉冲重复时间， 为等效收发阵元中心到目标的距离，,  为第1个等效收发阵元与Z轴夹角。则在第个快拍时刻，坐标系 中点坐标转变为坐标系 中的坐标。其中轴由和 的连线确定，轴和轴分别由右手定则确定； 轴由第个等效阵元和 的连线确定， 轴和 分别由右手定则确定。则经过平动补偿后，式(1)中的回波信号可以表示为sap(fb) =Q ∑q=1qexp[ j4πfbc(Rqap Rap)](2)Rap aOp其中， 表示第个等效阵元到转动中心 的距离。则Rqap Rap=√(Rap+ ^zqap)2+ (^yqap)2+ (^xqap)2 Rap=√(Rap)2+ 2Rap^zqap+(^rqap)2 Rap^zqap+(^rqap)22Rap(3)^rqap= (^xqap)2+ (^yqap)2+ (^zqap)2^rqap≪Rap Rqap Rap其中， ，在远场条件下近似认为 ，则根据文献[6]， 可以表示为Rqap Rap ^zqap=^xq0sin(φa) + ^ yq0cos(φa)sin(p)+ ^zq0cos(φa)cos(p) (4)将式(4)代入式(2)可得sap(fb)Q ∑q=1qexp[ j4πfcc(^xq0sin(φa)+ ^yq0cos(φa)sin(p)+^ zq0cos(φa)cos(p))](5) fcfb/cfc/c
0 0 d≪R0φa0 p0当发射信号带宽远小于中心频率时，可近似认为 ，假设 ，在远场条件下且观测时间较短时，则 ,  ,  。则式(5)可以近似为sap(b)Q ∑q=1qexp[ j4πfcc(^ zq0)]
exp[ j4π(b 1)∆fc(^ zq0)]
exp[ j4πfcc(^xq0(a 1)dR0+ ^yq0!(p 1)Tp)](6) apq= 4πfcc(^xq0(a 1)dR0+ ^yq0!(p 1)Tp)xx令 ，则通过对方向的相位求偏导数，可以得到沿方向相位变化为∆ x=d apqda=4πfcdcR0^xq0(7)x∆x假设两个散射中心在方向上的距离为 ，则总相位差可表示为 x=4πfcdcR0^xq0(MN 1) (8) x2π xx当满足 时，两个散射点在方向才能被分辨，所以方向上的分辨率为px= ∆xmin=cR02dfc(MN 1)(9) M1个发射阵元 N1个接收阵元等效收发节点 发射节点 接收节点2N1d d 2d 图 2 线性收发阵元的等效收发阵元示意图Fig. 2  The equivalent transceiver array element for linearreceiving and transmitting array elements YZXOqpjaO0ˆ xapˆ ˆ ˆ (xap,yap,zap)q q qˆ ˆ ˆ (x0,y0,z0)q q qˆ yapˆ zapqRapRapav0ˆ x0ˆ y0ˆ z0Op 图 3 目标移动示意图Fig. 3  Motion of the target第x期 邓理康等：一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法 3y同理可得方向的分辨率为py= ∆ymin=c2fc!PTp(10)zpz =c/(2Bw) Bwa;b;p而方向的分辨率为 ，其中 为信号带宽。所以对  3个维度分别做傅里叶变换就能得到三维图像。所以信号的张量模型可以表示为S=X1F12F23F3(11) n(n= 1;2;3)nS 2CAPBX 2CAPBF12CAAF22CPPF32CBB其中， ，表示张量和矩阵的模态积， 表示三维信号， 表示三维图像。 ,  ,  表示全傅里叶变换矩阵。当阵元数目减少，或由于噪声或硬件原因导致完整回波出现缺失时(相当于对回波的3个维度做稀疏采样)，如果继续采用傅里叶变换将会导致图像质量严重下降。因此我们利用图像的稀疏性，引入压缩感知方法对回波信号进行处理，用较少的数据就能恢复图像。但是传统的压缩感知算法需要将三维数据展开成一维形式。将式(11)展开成一维形式为s= (F3F2F1)x (12) s=vec(S); x= vec(X)606060 2160021600其中，符号表示克罗内克积，。假设三维信号的维度为 ，当转化为一维信号后，感知矩阵的维度为 ，这需要消耗巨大的内存和存储容量，限制了算法的进一步使用。所以针对上述问题，本文采用多维ADMM (MultiDimensional ADMM, MD-ADMM)算法对图像进行稀疏重构。3    多维ADMM稀疏成像方法3.1  多维ADMM算法下面以ADMM算法为基础，推导MD-ADMM算法，首先推导算法的三维张量形式，再推广到多维情形。当完整回波出现缺失时，式(11)可以表示为S=X1F(1)2F(2)3F(3)(13) S 2CMNKX 2CUVWUVW M; N;K F(1)2CMUF(2)2CNVF(3)2CKWF(1)=T1F1 F(2)=T2F2 F(3)=T3F3 F1;F2;F3 UU VV WW T1 2CMU;T2 2CNV; T3 2CKWG;H;J SG2[1;2;


;U]T; H2[1;2;


;V]T;J2[1;2;


;W]TM;N;K其中， 表示降采样信号，为待恢复图像。设完整信号维数为 , 分别为对信号3个维度的采样数。 ,,  ，分别为3个维度的部分傅里叶变换矩阵。令 ,  , ，其中 分别为维数为 ,  ,的全傅里叶矩阵。为傅里叶采样矩阵。令 分别表示对张量 3个维度的采样序列，其中，其序列长度分别为 。则T1(m;u) ={1; Gm=u0; 其他; T2(n;v) ={1; Hn=v0; 其他T3(k;w) ={1; Jk=w0; 其他  ,将式(13)展开成一维形式为s=(F(3)F(2)F(1))x (14)s=vec(S); x=vec(X) vec(
)sxx其中， ，其中 表示将张量向量化。在式(14)的信号模型中，直接从已知信号中重构出将有无穷多组解，属于病态问题，因此需要引入额外条件，压缩感知理论可以利用待恢复信号的稀疏性从无穷多组解中找出最稀疏的解。雷达图像通常由若干强散射点组成，因此待恢复信号的稀疏性一般情况下是成立的。理想情况下通常用l0范数表示信号的稀疏性，但是基于l0范数最小化的稀疏恢复问题不易求解，属于NP难问题。因此，一般采用l1范数替代l0范数，并且l1范数最小化问题通常可以转化为凸问题。本文采用ADMM算法[15]解决基于l1范数的最小化问题，该方法可以将高维的测量矩阵分解为张量的模态积，用张量元素除法替代了计算效率低的矩阵求逆，提高了计算效率。因此式(14)中，基于l1范数最小化优化问题可以表示为^ x=argminx{s (F(3)F(2)F(1))x22+∥x∥1}(15)z引入辅助变量，则原基于l1范数的最小化问题可以等价于以下带等式约束的优化问题^ x=argminx;z{s (F(3)F(2)F(1))x22+∥z∥1};s:t: x z= 0 (16)其增广拉格朗日函数可以写为L(x;z;
) =s (F(3)F(2)F(1))x22+∥z∥1+
H(x z) +2∥x z∥22(17)
其中，为对偶变量，为惩罚系数，ADMM算法将式(15)中的问题分解为如式(18)的3个子问题x(k+1)=argminxL(x;z(k);
(k))z(k+1)=argminzL(x(k+1);z;
(k))
(k+1)=
(k)+(x(k+1) z(k+1))9>=>;(18)4 雷  达  学  报 第10卷kkL(x;z;
)xzx;z
其中，上标表示第次迭代，式(18)中前两个子问题可以通过对函数 中的和分别求导并令导数为0求得，经求解式(18)中 以及的迭代式为x(k+1)=[(F(3)F(2)F(1))H
(F(3)F(2)F(1))+IUV W] 1
[(F(3)F(2)F(1))Hs+z(k)
(k)]z(k+1)=ST(x(k+1)+
(k);)
(k+1)=
(k)+(x(k+1) z(k+1))9>>=>>;(19) ST(x;a) = (x/jxj)max(jxj   x;0) F(2)=T2F2 F(3)=T3F3其中，ST为软阈值(Soft Threshold)函数，其表达式为 。将F(1)=T1·F1,  ,  代入式(19)可得x(k+1)=[(F3F2F1)H(B3B2B1)
(F3F2F)+IUV W] 1
[(F3F2F)H
(T3T2T1)Hs+z(k)
(k)](20)B1 =TH1T1;B2 =TH2T2;B3 =TH3T3其中， ，通过化简可得x(k+1)=(F3F2F1)H[(B3B2B1)+IUV W] 1[(T3T2T1)Hs+ (F3F2F1)(z(k)
(k))](21)将式(21)写成张量形式，可得X(k+1)={[(S1T12T23T3)+((Z(k)  A(k))1F12F23F3)]1FH12FH23FH3}⊘[G+1UVW](22) 1UVW UVW ⊘Gz(k+1)
(k+1)其中， 表示元素全为1、维数为的三维张量，表示张量的元素除法，表示对接收回波三维方向的采样，其值设为0或1，分别表示是否被采样到。同理式(19)中 和 的迭代式同样可以写成如式(23)和式(24)的张量形式Z(k+1)=ST(X(k+1)+
(k);)(23)A(k+1)=A(k)+(X(k+1)  Z(k+1))(24)XZA0联合迭代式(22)—式(24)，就可得到图像。初始参数设置如下：和的初值设定为，的值根据数据进行调整，的值取1。U1U2


Ui X(k+1)将三维形式进一步推广，假设多维张量的维数为 ，则张量 的更新表达式为X(k+1)={[(S1T12T2


iTi)+((Z(k)  A(k))1F12F2


iFi)]1FH12FH2


iFHi}⊘(G+1U1U2


Ui)Z(k+1)=ST(X(k+1)+
(k);)A(k+1)=A(k)+(X(k+1)  Z(k+1))9>>=>>;(25)3.2  计算复杂度分析XL2L3 ≫L1MD-ADMM算法中图像更新表达式(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]。在本实验中，迭代次数 ，且算法DR-2D-SL0的单次迭代计算复杂度最高，因此3种算法的计算复杂度排序为：MD-ADMM<MD-SL0<DR-2D-SL0。4    实验仿真4.1  仿真数据分析v0 = 200 R0= 10000 fc= 10 Bw= 150 B= 60 A= 60 PRF= 80 P= 60实验仿真数据设置如下：目标飞行速度 m/s，雷达距目标中心的距离  m，发射信号中心频率  GHz，带宽  MHz，发射步进频信号个数 。设收发阵列为10发6收MIMO线阵，等效收发阵元个数 ，等效收发阵元间距d=2.5 m。脉冲重复频率  Hz，快拍数 。仿真中使用的点散射模型如图4所示。H当回波数据完整时，对3个维度直接进行傅里叶变换后得到的图像如图5所示。由图5可知当回波数据完整时，直接对3个维度做傅里叶变换可以得到质量较高的三维图像。本文提取图5中的三维散点坐标，将散射点所在位置幅值设为1，图像中其余位置幅值设为0，形成参考三维图像。下面对第x期 邓理康等：一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法 5稀疏采样回波进行成像，为了获得三维稀疏回波，对回波进行稀疏采样，具体采样方式如图6所示。首先采用随机采样方式，对稀疏度(每个维度的稀疏度相同)分别为50%, 33.3%, 25%的回波进行三维成像处理，其中信号添加信噪比为20 dB的高斯白噪声。分别采用RD, MD-SL0, DR-2D-SL0,MD-ADMM 4种算法对目标进行三维成像(其中算法MD-SL0和算法DR-2D-SL0统称为SL0算法)，并采用三视图进行展示。在实验中本文调整参数让每一种算法都达到其最佳成像效果。图7—图9分别为稀疏度为50.0%, 33.3%, 25.0%的不同算法成像结果。由图7—图9可知，当回波稀疏时，由于受到旁瓣干扰，传统RD算法将会失效，得到的图像分辨率较低。由于利用了图像的稀疏性，采用了压缩感知方法，算法DR-2D-SL0, MD-SL0, MD-ADMM都得到了质量较高的图像。 20   10   0  -10  -20-20  -10   0   10   20y (m)x (m)(a) 目标三维点图(a) 3D point diagram of the target-15 -10 -5 0 5 10 15y (m)0246810z (m)0246810z (m)(c) 正视图(c) Front view-15 -10 -5 0 5 10 15-20-15-10-505101520y (m)x (m)(b) 俯视图(b) Top view-20 -10 0 10 200246810x (m)z (m)(d) 侧视图(d) Side view 图 4 仿真目标三维散点图Fig. 4  3D scatter of simulation target 俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元(a) 俯视图(a) Top view(b) 正视图(b) Front view方位单元(c) 侧视图(c) Side view 图 5 完整回波数据图像三视图Fig. 5  Three views of image with the complete echo6 雷  达  学  报 第10卷 (a) 随机采样(a) Random sampling(b) 块采样(b) Block sampling 图 6 回波采样形式Fig. 6  Undersampling masks of random sampling and block sampling 俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元(a) RD算法俯视图(a) Top view generated by RD algorithm(b) RD算法正视图(b) Front view generated by RD algorithm方位单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元方位单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元方位单元(c) RD算法侧视图(c) Side view generated by RD algorithm(d) SL0算法俯视图(d) Top view generated by SL0 algorithm(e) SL0算法正视图(e) Front view generated by SL0 algorithm(f) SL0算法侧视图(f) Side view generated by SL0 algorithm(g) MD-ADMM算法俯视图(g) Top view generated byMD-ADMM algorithm(h) MD-ADMM算法正视图(h) Front view generated byMD-ADMM algorithm(i) MD-ADMM算法侧视图(i) Side view generated byMD-ADMM algorithm 图 7 稀疏度为50%时图像Fig. 7  Image when sparsity is 50%第x期 邓理康等：一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法 7为了进一步定量比较4种算法，表1给出了随机稀疏采样条件下4种算法数值结果，其中包括图像熵、峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)以及计算时间。在图像处理领域，图像熵和PSNR可以一定程度上反映图像质量。在三维图像中图像熵的定义为式(26)E(X)= Sum{X2C⊙ln[X2C]}(26)Sum(
)C=∑Uu=1∑Vv=1∑Ww=1x2uvw U;V;WX其中， 表示张量的所有元素之和，,  分别为张量的XXH维数。假设对进行了归一化，则相对于参考图像的均方误差可以表示为MSE=1UV WU ∑u=1V ∑v=1W ∑w=1∥xuvw huvw∥ (27)huvw H其中， 为参考三维图像的元素。定义PSNR为PSNR= 10
lg(13/MSE)(28)由表1以及图7—图9可知，虽然RD算法计算效率最高，但是图像质量最差。由于参数设置一致，算法DR-2D-SL0和MD-SL0图像熵与PSNR相同， 俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060(a) RD算法俯视图(a) Top view generated by RD algorithm(b) RD算法正视图(b) Front view generated by RD algorithm方位单元距离单元10 20 30 40 50 60102030405060距离单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060方位单元距离单元10 20 30 40 50 60102030405060距离单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060方位单元距离单元10 20 30 40 50 60102030405060距离单元(c) RD算法侧视图(c) Side view generated by RD algorithm(d) SL0算法俯视图(d) Top view generated by SL0 algorithm(e) SL0算法正视图(e) Front view generated by SL0 algorithm(f) SL0算法侧视图(f) Side view generated by SL0 algorithm(g) MD-ADMM算法正视图(g) Top view generated byMD-ADMM algorithm(h) MD-ADMM算法侧视图(h) Front view generated byMD-ADMM algorithm(i) MD-ADMM算法侧视图(i) Side view generated byMD-ADMM algorithm 图 8 稀疏度为33.3%时的图像Fig. 8  Image when sparsity is 33.3%8 雷  达  学  报 第10卷但是MD-SL0算法直接对三维张量进行处理，而DR-2D-SL0要将三维张量展开成二维矩阵，这增加了计算量，所以MD-SL0的计算时间更短，计算效率更高。与其余基于压缩感知算法相比，所提MD-ADMM算法在不同稀疏度下的图像熵最小，峰值信噪比最大，并且计算时间最短，这验证了所提算法的有效性。接着采用块稀疏采样方式，对稀疏度(每个维度的稀疏度相同)为50%, 33.3%, 25%的回波进行成像处理。图10为采用不同算法对信噪比为20 dB，稀疏度为50%的回波进行处理后得到的图像。表2为块稀疏采样条件下不同算法的数值结果，由图10和表2可知，所提MD-ADMM算法在不同稀疏度的块稀疏采样条件下得到的图像熵最小，峰值信噪比最大，并且计算时间最短，这验证了所提算法的有效性。下面比较不同算法在相同稀疏度(3个维度的稀疏度相同)，不同信噪比下的成像性能，其中回波信号采用随机采样方式，每个维度的稀疏度均为25%，并添加均值为0的高斯白噪声。图11—图13分别为在信噪比为–5, 0, 10 dB条件下不同算法得到的目标三视图。表3为不同信噪比条件下的数值 俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元(a) RD算法俯视图(a) Top view generated by RD algorithm(b) RD算法正视图(b) Front view generated by RD algorithm方位单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元方位单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元方位单元(c) RD算法侧视图(c) Side view generated by RD algorithm(d) SL0算法俯视图(d) Top view generated by SL0 algorithm(e) SL0算法正视图(e) Front view generated by SL0 algorithm(f) SL0算法侧视图(f) Side view generated by SL0 algorithm(g) MD-ADMM算法正视图(g) Top view generated byMD-ADMM algorithm(h) MD-ADMM算法侧视图(h) Front view generated byMD-ADMM algorithm(i) MD-ADMM算法侧视图(i) Side view generated byMD-ADMM algorithm 图 9 稀疏度为25%时的图像Fig. 9  Image when sparsity is 25%第x期 邓理康等：一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法 9表 1 随机稀疏采样条件下数值结果Tab. 1   Numerical results under random sparse sampling condition稀疏度 算法 图像熵 PSNR 计算时间50.0%RD 9.490 34.075 0.012DR-2D-SL0 3.138 57.714 35.053MD-SL0 3.138 57.714 6.954MD-ADMM 3.006 61.678 3.51833.3%RD 10.505 28.430 0.010DR-2D-SL0 3.200 52.702 29.091MD-SL0 3.200 52.702 12.546MD-ADMM 2.959 53.888 6.08825.0%RD 10.965 26.003 0.0131DR-2D-SL0 3.373 48.755 54.008MD-SL0 3.373 48.755 28.509MD-ADMM 3.023 49.386 11.996 俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元(a) RD算法俯视图(a) Top view generated by RD algorithm(b) RD算法正视图(b) Front view generated by RD algorithm方位单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元方位单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元方位单元(c) RD算法侧视图(c) Side view generated by RD algorithm(d) SL0算法俯视图(d) Top view generated by SL0 algorithm(e) SL0算法正视图(e) Front view generated by SL0 algorithm(f) SL0算法侧视图(f) Side view generated by SL0 algorithm(g) MD-ADMM算法俯视图(g) Top view generated byMD-ADMM algorithm(h) MD-ADMM算法正视图(h) Front view generated byMD-ADMM algorithm(i) MD-ADMM算法侧视图(i) Side view generated byMD-ADMM algorithm 图 10 稀疏度为50%的块稀疏采样图像Fig. 10  The image of a sparsity of 50% by block sampling10 雷  达  学  报 第10卷表 2 块稀疏采样条件下数值结果Tab. 2   Numerical results under block sparse sampling condition稀疏度 算法 图像熵 PSNR 计算时间50.0%RD 8.336 32.885 0.012DR-2D-SL0 3.246 51.232 41.907MD-SL0 3.246 51.232 8.029MD-ADMM 3.037 53.274 3.66033.3%RD 8.553 27.495 0.0098DR-2D-SL0 3.700 42.188 42.518MD-SL0 3.700 42.188 17.954MD-ADMM 3.233 42.336 5.75825.0%RD 9.500 26.906 0.0092DR-2D-SL0 3.894 45.196 51.463MD-SL0 3.894 45.196 25.702MD-ADMM 3.337 47.201 11.536 俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元(a) RD算法俯视图(a) Top view generated by RD algorithm(b) RD算法正视图(b) Front view generated by RD algorithm方位单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元方位单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元方位单元(c) RD算法侧视图(c) Side view generated by RD algorithm(d) SL0算法俯视图(d) Top view generated by SL0 algorithm(e) SL0算法正视图(e) Front view generated by SL0 algorithm(f) SL0算法侧视图(f) Side view generated by SL0 algorithm(g) MD-ADMM算法俯视图(g) Top view generated byMD-ADMM algorithm(h) MD-ADMM算法正视图(h) Front view generated byMD-ADMM algorithm(i) MD-ADMM算法侧视图(i) Side view generated byMD-ADMM algorithm 图 11 稀疏度为25%信噪比为–5 dB的图像Fig. 11  The image when sparsity is 25% and SNR=–5 dB第x期 邓理康等：一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法 11结果。由图11—图13和表3可知，在稀疏孔径条件下RD算法基本失效。在不同信噪比条件下，所提MD-ADMM算法图像熵最小，PSNR最大，并且计算时间最短，这证明了所提算法对噪声的鲁棒性最强。4.2  实测数据分析MIMO雷达实验系统目前还在搭建中，还存在发射信号带宽过窄、收发阵列同步性差等问题，导致MIMO-ISAR回波目前还难以获取，也是下一步要着重解决的问题。因此采用二维Yak-42飞机IS-AR实测数据，以验证所提MD-ADMM算法在有限采样数据条件下的有效性。实验数据参数设置如下：雷达发射信号的载频为5520 MHz，信号带宽为400 MHz，快时间采样频率为10 MHz，脉冲宽度为25.6 µs，观测目标为Yak-42飞机。假设接收信号已经做了包络对齐、平动补偿以及自聚焦，共接收到256个脉冲，每个脉冲信号包含256个快时间采样。本文采样稀疏采样方式，抽取96个脉冲，以及128个快时间信号。图14为不同算法对二维稀疏信号成像结果。图15为图14信号中添加0 dB的高斯白噪声的结果。表4为不同算法在不同信噪比条件下的数值结果。由图14、图15以及表4可知，RD算 俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元(a) RD算法俯视图(a) Top view generated by RD algorithm(b) RD算法正视图(b) Front view generated by RD algorithm方位单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元方位单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元方位单元(c) RD算法侧视图(c) Side view generated by RD algorithm(d) SL0算法俯视图(d) Top view generated by SL0 algorithm(e) SL0算法正视图(e) Front view generated by SL0 algorithm(f) SL0算法侧视图(f) Side view generated by SL0 algorithm(g) MD-ADMM算法俯视图(g) Top view generated byMD-ADMM algorithm(h) MD-ADMM算法正视图(h) Front view generated byMD-ADMM algorithm(i) MD-ADMM算法侧视图(i) Side view generated byMD-ADMM algorithm 图 12 稀疏度为25%信噪比为0 dB图像Fig. 12  The image when sparsity is 25% and SNR=0 dB12 雷  达  学  报 第10卷表 3 不同信噪比条件下数值结果Tab. 3   Numerical results under different SNR conditions信噪比 算法 图像熵 PSNR 计算时间–5 dBRD 11.8042 18.320 0.0177DR-2D-SL0 7.3218 37.639 49.998MD-SL0 7.3218 37.639 23.540MD-ADMM 4.779 43.352 10.7750 dBRD 11.535 23.452 0.013DR-2D-SL0 5.4086 43.285 48.776MD-SL0 5.4086 43.285 23.303MD-ADMM 2.977 47.449 10.71510 dBRD 11.116 27.145 0.0124DR-2D-SL0 3.948 45.619 48.754MD-SL0 3.948 45.619 23.484MD-ADMM 3.066 47.561 10.521 俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元(a) RD算法俯视图(a) Top view generated by RD algorithm(b) RD算法正视图(b) Front view generated by RD algorithm方位单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元方位单元俯仰单元方位单元10 20 30 40 50 60102030405060俯仰单元10 20 30 40 50 60102030405060距离单元10 20 30 40 50 60102030405060距离单元方位单元(c) RD算法侧视图(c) Side view generated by RD algorithm(d) SL0算法俯视图(d) Top view generated by SL0 algorithm(e) SL0算法正视图(e) Front view generated by SL0 algorithm(f) SL0算法侧视图(f) Side view generated by SL0 algorithm(g) MD-ADMM算法俯视图(g) Top view generated byMD-ADMM algorithm(h) MD-ADMM算法正视图(h) Front view generated byMD-ADMM algorithm(i) MD-ADMM算法侧视图(i) Side view generated byMD-ADMM algorithm 图 13 稀疏度为25%信噪比为10 dB的图像Fig. 13  The image when sparsity is 25% and SNR=10 dB第x期 邓理康等：一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法 13 多普勒单元距离单元50 100 150 200 250501001502002500.10.20.30.40.50.60.70.80.9多普勒单元距离单元500000 100 150 200 250501001502002500.10.20.30.40.50.60.70.80.9(a) 完整数据二维FFT结果(a) Two dimensional FFT result of full data(c) MD-SL0成像结果(c) Imaging result via MD-SL0多普勒单元距离单元50 100 150 200 250501001502002500.10.20.30.40.50.60.70.80.9多普勒单元距离单元50 100 150 200 250501001502002500.10.20.30.40.50.60.70.80.9(b) RD成像结果(b) Imaging result via RD(d) MD-ADMM的成像结果(d) Imaging result via MD-ADMM 图 14 实测数据结果Fig. 14  Measured ISAR data results 距离单元5050 100 150 200 2501001502002500.10.20.30.40.50.60.70.80.9(a) 完整数据二维FFT结果(a) Two dimensional FFT result of full data(c) MD-SL0成像结果(c) Imaging result via MD-SL0多普勒单元距离单元50500000 100 150 200 2501001502002500.10.20.30.40.50.60.70.80.9多普勒单元(b) RD成像结果(b) Imaging result via RD(d) MD-ADMM的成像结果(d) Imaging result via MD-ADMM距离单元5050 100 150 200 2501001502002500.10.20.30.40.50.60.70.80.9多普勒单元距离单元5050 100 150 200 2501001502002500.10.20.30.40.50.60.70.80.9多普勒单元 图 15 信噪比为0 dB实测数据结果Fig. 15  Measured ISAR data results under SNR=0 dB14 雷  达  学  报 第10卷法基本失效，算法MD-SL0, MD-ADMM在二维稀疏采样条件下依然能够得到清晰图像，但相比MD-SL0算法，MD-ADMM算法得到的图像熵最小，峰值信噪比最大，并且计算时间最短，进一步验证了所提算法的有效性。5    结论本文提出一种多维ADMM稀疏恢复算法，本算法可以用于恢复多维稀疏信号，而无需将多维信号转换为一维，这极大地减少了存储量和计算量。通过该算法，实现了一种实用的MIMO-ISAR成像。与其他算法相比，本算法具有存储容量小、计算效率高、成像质量好的优点。仿真和实测数据结果均证明了本算法的有效性。参 考 文 献MUSMAN S, KERR D, and BACHMANN C. Automaticrecognition of ISAR ship images[J]. IEEE Transactions onAerospace and Electronic Systems, 1996, 32(4): 1392–1404.doi: 10.1109/7.543860.[1]MARTORELLA M, GIUSTI E, DEMI L, et al. Targetrecognition by means of polarimetric ISAR images[J]. IEEETransactions on Aerospace and Electronic Systems, 2011,47(1): 225–239. doi: 10.1109/TAES.2011.5705672.[2]BARANIUK R and STEEGHS P. Compressive radarimaging[C]. 2007 IEEE Radar Conference, Waltham, USA,2007: 128–133.[3]DUAN Guangqing, WANG Dangwei, MA Xiaoyan, et al.Three-dimensional imaging via wideband MIMO radarsystem[J]. IEEE Geoscience and Remote Sensing Letters,2010, 7(3): 445–449. doi: 10.1109/LGRS.2009.2038728.[4]HU Xiaowei, TONG Ningning, WANG Heming, et al. [5]Multiple-input-multiple-output radar superresolution three-dimensional imaging based on multidimensional smoothedL0[J]. Journal of Applied Remote Sensing, 2016, 10(3):035017. doi: 10.1117/1.JRS.10.035017.WANG Yong and LI Xuelu. 3-D imaging based oncombination of the ISAR technique and a MIMO radarsystem[J]. IEEE Transactions on Geoscience and RemoteSensing, 2018, 56(10): 6033–6054. doi: 0.1109/TGRS.2018.2829912.[6]CANDES E J and WAKIN M B. An introduction tocompressive sampling[J]. IEEE Signal Processing Magazine,2008, 25(2): 21–30. doi: 10.1109/MSP.2007.914731.[7]DONOHO D L. Compressed sensing[J]. IEEE Transactionson Information Theory, 2006, 52(4): 1289–1306. doi: 10.1109/TIT.2006.871582.[8]DUARTE M F and BARANIUK R G. Kroneckercompressive sensing[J]. IEEE Transactions on ImageProcessing, 2012, 21(2): 494–504. doi: 10.1109/TIP.2011.2165289.[9]QIU Wei, MARTORELLA M, ZHOU Jianxiong, et al.Three-dimensional inverse synthetic aperture radar imagingbased on compressive sensing[J]. IET Radar, Sonar &Navigation, 2015, 9(4): 411–420.[10]HU Xiaowei, TONG Ningning, GUO Yiduo, et al. MIMOradar 3-D imaging based on multi-dimensional sparserecovery and signal support prior information[J]. IEEESensors Journal, 2018, 18(8): 3152–3162. doi: 10.1109/JSEN.2018.2810705.[11]QIU Wei, ZHOU Jianxiong, ZHAO Hongzhong, et al.Three-dimensional sparse turntable microwave imagingbased on compressive sensing[J]. IEEE Geoscience andRemote Sensing Letters, 2015, 12(4): 826–830. doi:10.1109/LGRS.2014.2363238.[12]ZHANG Shuanghui, LIU Yongxiang, and LI Xiang.Computationally  efficient  sparse  aperture  ISARautofocusing and imaging based on fast ADMM[J]. IEEETransactions on Geoscience and Remote Sensing, 2020,58(12): 8751–8765. doi: 10.1109/TGRS.2020.2990445.[13]ZHU Yutao, SU Yi, and YU Wenxian. An ISAR imagingmethod based on MIMO technique[J]. IEEE Transactionson Geoscience and Remote Sensing, 2010, 48(8): 3290–3299.doi: 10.1109/TGRS.2010.2045230.[14]BOYD S, PARIKH N, CHU E, et al. Distributedoptimization and statistical learning via the alternatingdirection method of multipliers[J]. Foundations andTrends® in Machine Learning, 2011, 3(1): 1–122.[15]表 4 实测数据不同信噪比条件下的数值结果Tab. 4  Numerical results of measured data for differentsignal-to-noise ratio conditions信噪比 算法 图像熵 PSNR 计算时间—RD 9.639 29.816 0.010MD-SL0 6.459 40.378 10.553MD-ADMM 5.296 41.364 3.7430 dBRD 10.242 26.909 0.0131MD-SL0 7.571 35.241 12.064MD-ADMM 6.183 39.090 4.274第x期 邓理康等：一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法 15作 者 简 介 邓理康(1991–)，男，福建建阳人，国防科技大学电子科学学院在读研究生，研究方向为双站雷达成像、MIMO雷达成像。E-mail: denglikang09@126.com张双辉(1989–)，男，湖南长沙人，博士，国防科技大学电子科学学院副研究员，研究方向为雷达成像、压缩感知、贝叶斯推断。E-mail: shzhang3@126.com张　弛(1994–)，男，湖北孝感人，国防科技大学电子科学学院在读博士生，研究方向为雷达成像、压缩感知、贝叶斯学习。E-mail: zcnudt@126.com刘永祥(1976–)，男，河北唐山人，博士，国防科技大学电子科学学院教授，博士生导师，研究方向为目标微动特性分析与识别。E-mail: lyx_bible@sina.com16 雷  达  学  报 第10卷