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

邓理康 张双辉 张弛 刘永祥

邓理康, 张双辉, 张弛, 等. 一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法[J]. 雷达学报, 待出版. doi:  10.12000/JR20132
引用本文: 邓理康, 张双辉, 张弛, 等. 一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法[J]. 雷达学报, 待出版. doi:  10.12000/JR20132
DENG Likang, ZHANG Shuanghui, ZHANG Chi, et al. A multiple-input multiple-output inverse synthetic aperture radar imaging method based on multidimensional alternating direction method of multipliers[J]. Journal of Radars, in press. doi:  10.12000/JR20132
Citation: DENG Likang, ZHANG Shuanghui, ZHANG Chi, et al. A multiple-input multiple-output inverse synthetic aperture radar imaging method based on multidimensional alternating direction method of multipliers[J]. Journal of Radars, in press. doi:  10.12000/JR20132

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

doi: 10.12000/JR20132
基金项目: 国家自然科学基金(61801484, 61921001)
详细信息
    作者简介:

    邓理康(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.com

    通讯作者:

    张双辉 shzhang3@126.com

  • 责任主编:张弓 Corresponding Editor: ZHANG Gong
  • 中图分类号: TN957.51

A Multiple-Input Multiple-Output Inverse Synthetic Aperture Radar Imaging Method Based on Multidimensional Alternating Direction Method of Multipliers

Funds: The National Natural Science Foundation of China (61801484, 61921001)
More Information
图(15) / 表 (4)
计量
  • 文章访问数:  55
  • HTML全文浏览量:  9
  • PDF下载量:  15
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-10-19
  • 修回日期:  2021-01-27
  • 网络出版日期:  2021-02-08

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

doi: 10.12000/JR20132
    基金项目:  国家自然科学基金(61801484, 61921001)
    作者简介:

    邓理康(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.com

    通讯作者: 张双辉 shzhang3@126.com
  • 责任主编:张弓 Corresponding Editor: ZHANG Gong
  • 中图分类号: TN957.51

摘要: 基于傅里叶变换的传统逆合成孔径雷达(ISAR)成像方法存在数据存储量大、数据采集时间长的问题。压缩感知(CS)理论利用图像的稀疏性,可以利用有限的数据恢复图像,这极大降低了数据采集成本。但对于多维数据,传统压缩感知方法要将多维数据转化成一维向量,这造成了很大存储和计算负担。因此,该文提出一种基于多维度-交替方向乘子法(MD-ADMM)的多输入多输出-逆合成孔径雷达(MIMO-ISAR)成像快速稀疏重建方法。首先建立基于张量信号的压缩感知模型,然后用ADMM算法对模型进行优化,将测量矩阵分解为张量模态积,用张量元素除法替代矩阵求逆,显著减少所需的内存和计算负担。该方法只需少量的数据采样,就能实现快速成像。与其他基于张量的压缩感知方法相比,该方法具有鲁棒性强、图像质量好、计算效率高的优点。仿真和实测数据验证了该方法的有效性。

注释:
1)  责任主编:张弓 Corresponding Editor: ZHANG Gong

English Abstract

邓理康, 张双辉, 张弛, 等. 一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法[J]. 雷达学报, 待出版. doi:  10.12000/JR20132
引用本文: 邓理康, 张双辉, 张弛, 等. 一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法[J]. 雷达学报, 待出版. doi:  10.12000/JR20132
DENG Likang, ZHANG Shuanghui, ZHANG Chi, et al. A multiple-input multiple-output inverse synthetic aperture radar imaging method based on multidimensional alternating direction method of multipliers[J]. Journal of Radars, in press. doi:  10.12000/JR20132
Citation: DENG Likang, ZHANG Shuanghui, ZHANG Chi, et al. A multiple-input multiple-output inverse synthetic aperture radar imaging method based on multidimensional alternating direction method of multipliers[J]. Journal of Radars, in press. doi:  10.12000/JR20132
    • 逆合成孔径雷达(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 Reduction 2D Smoothed l0-norm, DR-2D-SL0)重构方法求解稀疏优化问题,得到散射体的二维图像。最后,通过将得到的二维图像重新还原成三维张量来得到三维图像。文献[11,12]首先将信号表示为张量和矩阵的模态积形式,然后采用一种稀疏多维度平滑l0算法(MultiDimensional Smoothed l0-norm, MD-SL0)恢复三维图像,由于该算法能够直接处理三维图像,所以具有较低的计算复杂度和内存消耗,能够有效地进行三维重建。但是该方法需要调整的参数较多,并且在低信噪比条件下图像恢复效果一般。文献[13]将ADMM算法融入稀疏孔径成像ISAR自聚焦中,利用部分傅里叶矩阵字典将矩阵求逆转化为元素除法,并将二维ISAR图像的距离单元更新变为整体更新,在数秒内就能获得聚焦良好的图像。

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

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

      图  1  成像场景图

      Figure 1.  Geometry of imaging

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

      图  2  线性收发阵元的等效收发阵元示意图

      Figure 2.  The equivalent transceiver array element for linear receiving and transmitting array elements

      假设发射阵列发射步进频信号,发射频率${f_b} = {f_{c}} + \left( {b - 1} \right)\Delta f$,其中${f_c}$为发射信号中心频率,$\Delta f$为频率步长。经过信号分选后,第$a$个等效收发阵元在第$p$个快拍时刻的接收信号可以表示为

      $${s_a}_p({f_b}) = \sum\limits_{q = 1}^Q {{\sigma _q}} {\kern 1pt} \exp \left[ { - {\rm{j}}\frac{{4\pi {f_b}}}{\rm{c}} R_{ap}^q} \right]{\kern 1pt} $$ (1)

      其中,$a = 1,2, \cdots , A$表示等效收发阵元序号,A为等效阵元总数;$q = 1,2, \cdots , Q$表示散射点序号,$Q$为散射点总数;$p = 1,2, \cdots ,P$表示快拍序号,$P$为快拍数;${\rm{c}}$为光速;${\sigma _q}$为第$q$个散射点的散射强度;$R_{ap}^q$表示第$a$个等效收发阵元在第p个脉冲时刻到第q个散射点的距离。

      图3所示,目标在ZOY 面沿着Y 轴方向以速度${v_0}$作匀速直线运动,假设初始快拍时刻,目标参考点位于${O_0}$,第p个快拍时刻参考点位于${O_p}$, ${\varphi _a} \approx (a - 1)d/{R_0}\;{\rm{ + }}\;{\alpha _0}$为第$a$个等效阵元与Z 轴夹角,${\theta _p} \approx \omega (p - 1){T_p}$$O{O_p}$Z 轴夹角。其中$d$为等效阵元间的距离,$\omega $为目标运动等效转速,${T_p}$为脉冲重复时间,${R_0}$为等效收发阵元中心到目标的距离,$\omega \approx v/{R_0}$, ${\alpha _0}$为第1个等效收发阵元与Z轴夹角。则在第$p$个快拍时刻,坐标系${O_0} {\text{-}}{\hat x_0}{\hat y_0}{\hat z_0}$$q$点坐标$(\hat x_0^q,\hat y_0^q,\hat z_0^q)$转变为坐标系${O_p} {\text{-}} {\hat x_{ap}}{\hat y_{ap}}{\hat z_{ap}}$中的坐标$(\hat x_{ap}^q,\hat y_{ap}^q,\hat z_{ap}^q)$。其中${\hat z_0}$轴由$O$${O_0}$的连线确定,${\hat x_0}$轴和${\hat y_0}$轴分别由右手定则确定;${\hat z_{ap}}$轴由第$a$个等效阵元和${O_p}$的连线确定,${\hat x_{ap}}$轴和${\hat y_{ap}}$分别由右手定则确定。则经过平动补偿后,式(1)中的回波信号可以表示为

      图  3  目标移动示意图

      Figure 3.  Motion of the target

      $${s_a}_p({f_b}) = \sum\limits_{q = 1}^Q {{\sigma _q}} {\kern 1pt} \exp \left[ { - {\rm{j}}\frac{{4\pi {f_b}}}{{\rm{c}}} \left( {R_{ap}^q - {R_{ap}}} \right)} \right]$$ (2)

      其中,${R_{ap}}$表示第$a$个等效阵元到转动中心${O_p}$的距离。则

      $$ \begin{split} R_{ap}^q - {R_{ap}} &= \sqrt {{{\left( {{R_{ap}} + \hat z_{ap}^q} \right)}^2} + {{\left( {\hat y_{ap}^q} \right)}^2} + {{\left( {\hat x_{ap}^q} \right)}^2}} - {R_{ap}} \\ &= \sqrt {{{\left( {{R_{ap}}} \right)}^2} + 2{R_{ap}}\hat z_{ap}^q{\rm{ + }}{{\left( {\hat r_{ap}^q} \right)}^2}} - {R_{ap}} \\ &{\kern 1pt} \approx \hat z_{ap}^q + \frac{{{{\left( {\hat r_{ap}^q} \right)}^2}}}{{2{R_{ap}}}}\\[-15pt] \end{split} $$ (3)

      其中,$\hat r_{ap}^q = {(\hat x_{ap}^q)^2} + {(\hat y_{ap}^q)^2} + {(\hat z_{ap}^q)^2}$,在远场条件下近似认为$\hat r_{ap}^q \ll {R_{ap}}$,则根据文献[6],$R_{ap}^q - {R_{ap}}$可以表示为

      $$\begin{split} R_{ap}^q - {R_{ap}}{\kern 1pt} \approx &\hat z_{ap}^q \\ = & \hat x_0^q\sin({\varphi _a}) + \hat y_0^q\cos({\varphi _a})\sin({\theta _p}) \\ & + \hat z_0^q\cos({\varphi _a})\cos({\theta _p}) \end{split} $$ (4)

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

      $$\begin{split} {s_a}_p({f_b}) \approx &\sum\limits_{q = 1}^Q {{\sigma _q}} {\kern 1pt} \exp\left[ { - {\rm{j}}\frac{{4\pi {f_{\rm{c}}}}}{{\rm{c}}} \Big( {\hat x_0^q\sin({\varphi _a})} } \right. \\ & + \hat y_0^q\cos({\varphi _a})\sin({\theta _p}). { { + \hat z_0^q\cos({\varphi _a})\cos({\theta _p})} \Big)} ] \end{split} $$ (5)

      当发射信号带宽远小于中心频率${f_{\rm{c}}}$时,可近似认为${f_b}/{\rm{c}} \approx {f_{\rm{c}}}/{\rm{c}}$,假设${\alpha _0} \approx 0$,在远场条件下且观测时间较短时,则$d \ll {R_0}$, ${\varphi _a} \approx 0$, ${\theta _p} \approx 0$。则式(5)可以近似为

      $$ {\begin{split} {s_a}_p(b) \approx & \sum\limits_{q = 1}^Q {{\sigma _q}} {\kern 1pt} \exp\left[ { - {\rm{j}}\frac{{4\pi {f_c}}}{{\rm{c}}}\left( {\hat z_0^q} \right)} \right] \\ &\cdot \exp\left[ { - {\rm{j}}\frac{{4\pi (b - 1)\Delta f}}{{\rm{c}}}\left( {\hat z_0^q} \right)} \right] \\ &\cdot \exp\left[ { - {\rm{j}}\frac{{4\pi {f_c}}}{{\rm{c}} }\left( {\hat x_0^q\frac{{(a - 1)d}}{{{R_0}}} + \hat y_0^q\omega (p - 1){T_p}} \right)} \!\!\right] \end{split}} $$ (6)

      ${\psi _{apq}} = - \dfrac{{4\pi {f_{\rm{c}}}}}{{\rm{c}}}\! \!\left(\! {\hat x_0^q\dfrac{{(a - 1)d}}{{{R_0}}} + \hat y_0^q\omega (p - 1){T_p}} \!\right)$,则通过对$x$方向的相位求偏导数,可以得到沿$x$方向相位变化为

      $$\Delta {\psi _x} = \frac{{{\rm{d}}{\psi _{apq}}}}{{{\rm{d}}a}}{\rm{ = }}\frac{{4\pi {f_{\rm{c}}}d}}{{{\rm{c}}{R_0}}}\hat x_0^q$$ (7)

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

      $${\psi _x} = \frac{{4\pi {f_{\rm{c}}}d}}{{{\rm{c}}{R_0}}}\hat x_0^q\left( {MN - 1} \right)$$ (8)

      当满足${\psi _x} \ge 2\pi$时,两个散射点在$x$方向才能被分辨,所以$x$方向上的分辨率为

      $${p_x} = \Delta {x_{\min }} = \frac{{c{R_0}}}{{2d{f_{\rm{c}}}\left( {MN - 1} \right)}}$$ (9)

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

      $${p_y} = \Delta {y_{\min }} = \frac{{\rm{c}}}{{2{f_{\rm{c}}}\omega P{T_p}}}$$ (10)

      $z$方向的分辨率为${p_z} = {\rm{c}}/(2{B_{\rm{w}}})$,其中${B_{\rm{w}}}$为信号带宽。所以对$a,b,p$ 3个维度分别做傅里叶变换就能得到三维图像。所以信号的张量模型可以表示为

      $${\cal{S}} = {\cal{X}}{ \times _1}{{{F}}_1}{ \times _2}{{{F}}_2}{ \times _3}{{{F}}_3}$$ (11)

      其中,${ \times _n}\left( {n = 1,2,3} \right)$,表示张量和矩阵的$n$模态积,${\cal{S}} \in {{\mathbb{C}}^{A \times P \times B}}$表示三维信号,${\cal{X}} \in {{\mathbb{C}}^{A \times P \times B}}$表示三维图像。${{{F}}_1}\! \in \!{{\mathbb{C}}^{A \times A}}$, ${{{F}}_2} \!\in \!{{\mathbb{C}}^{P \times P}}$, ${{{F}}_3}\! \in \!{{\mathbb{C}}^{B \times B}}$表示全傅里叶变换矩阵。

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

      $$ {{s}} = \left( {{{{F}}_3} \otimes {{{F}}_2} \otimes {{{F}}_1}} \right){{x}} $$ (12)

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

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

      $${\cal{S}} = {\cal{X}}{ \times _1}{{{F}}^{\left( 1 \right)}}{ \times _2}{{{F}}^{\left( 2 \right)}}{ \times _3}{{{F}}^{(3)}}$$ (13)

      其中,${\cal{S}}\! \in \!{{\mathbb{C}}^{M \times N \times K}}$表示降采样信号,${\cal{X}}\! \in \!{{\mathbb{C}}^{U \times V \times W}}$为待恢复图像。设完整信号维数为$U \times V \times W$, $M, N,K$分别为对信号3个维度的采样数。${{{F}}^{\left( 1 \right)}} \in {{\mathbb{C}}^{M \times U}}$,${{{F}}^{\left( 2 \right)}} \in {{\mathbb{C}}^{N \times V}}$, ${{{F}}^{\left( 3 \right)}} \in {{\mathbb{C}}^{K \times W}}$,分别为3个维度的部分傅里叶变换矩阵。令${{{F}}^{(1)}}\! = \!{{{T}}_1}{{{F}}_1}$, ${{{F}}^{(2)}} = {{{T}}_2}{{{F}}_2}$, ${{{F}}^{(3)}} = {{{T}}_3}{{{F}}_3}$,其中${{{F}}_1},{{{F}}_2},{{{F}}_3}$分别为维数为$U \times U$, $V \times V$, $W \times W$的全傅里叶矩阵。${{{T}}_1} \in {{\mathbb{C}}^{M \times U}},{{{T}}_2} \in {{\mathbb{C}}^{N \times V}},$ ${{{T}}_3} \in {{\mathbb{C}}^{K \times W}}$为傅里叶采样矩阵。令${{G}},{{H}},{{J}}$分别表示对张量${\cal{S}}$ 3个维度的采样序列,其中,${{G}} \in {\left[ {1,2, \cdots,U} \right]^{\rm{T}}}, {{H}} \in {\left[ {1,2, \cdots,V} \right]^{\rm{T}}},{{J}} \in {\left[ {1,2, \cdots,W} \right]^{\rm{T}}}$其序列长度分别为$M,N,K$。则$ {T}_{1}\left(m,u\right)=\left\{\!\! \begin{array}{l}1,\;{{G}}_{m}=u\\ 0,\;{\text{其他}}\end{array} \right., \; $ ${T}_{2}\left(n,v\right)= \left\{ \begin{array}{l}1,\;\;{{H}}_{n}=v\\ 0,\;\;{\text{其他}}\end{array} \right. $, $ {T}_{3}\left(k,w\right)=\left\{ \begin{array}{l}1,\;\;{{J}}_{k}=w\\ 0,\;\;{\text{其他}}\end{array} \right.$将式(13)展开成一维形式为

      $${{s}} = \left({{{F}}^{\left( 3 \right)}} \otimes {{{F}}^{(2)}} \otimes {{{F}}^{\left( 1 \right)}}\right){{x}}$$ (14)

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

      $$\begin{split} \hat {{x}} = &\arg\mathop {\min}\limits_{{x}} \left\{ {\left\| {{{s}} - \left({{{F}}^{\left( 3 \right)}} \otimes {{{F}}^{(2)}} \otimes {{{F}}^{\left( 1 \right)}}\right){{x}}} \right\|_2^2} \right. \\ & . { + \lambda {{\left\| {{x}} \right\|}_1}} \} \end{split} $$ (15)

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

      $$\begin{split} &\hat {{x}} = \arg \mathop {\min }\limits_{{{x}},{{z}}} \left\{ {\left\| {{{s}} - \left({{{F}}^{\left( 3 \right)}} \otimes {{{F}}^{(2)}} \otimes {{{F}}^{\left( 1 \right)}}\right){{x}}} \right\|_2^2} \right. \\ &\quad\quad. { + \lambda {{\left\| {{z}} \right\|}_1} } \}, {\rm{s}}{\rm{.t}}{\rm{. }} \;\;\;{{x}} - {{z}} = 0 \end{split} $$ (16)

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

      $$\begin{split} {L_\rho }({{x}},{{z}},{{\alpha}} ) = &\left\| {{{s}} - \left({{{F}}^{\left( 3 \right)}} \otimes {{{F}}^{(2)}} \otimes {{{F}}^{\left( 1 \right)}}\right){{x}}} \right\|_2^2 \\ &+ \lambda {\left\| {{z}} \right\|_1} + {{{\alpha}} ^{\rm{H}}}({{x}} - {{z}}) + \frac{\rho }{2}\left\| {{{x}} - {{z}}} \right\|_2^2 \end{split} $$ (17)

      其中,${{\alpha}} $为对偶变量,$\rho $为惩罚系数,ADMM算法将式(15)中的问题分解为如式(18)的3个子问题

      $$\left. \begin{aligned} &{{{x}}^{\left( {k + 1} \right)}} = \arg \mathop {\min }\limits_{{x}} {L_\rho }\left({{x}},{{{z}}^{(k)}},{{{\alpha}} ^{(k)}}\right) \\ &{{{z}}^{(k + 1)}} = \arg \mathop {\min }\limits_{{z}} {L_\rho }\left({{{x}}^{(k + 1)}},{{z}},{{{\alpha}} ^{(k)}}\right) \\ &{{{\alpha}} ^{(k + 1)}} = {{{\alpha}} ^{(k)}} + \rho \left( {{{{x}}^{\left( {k + 1} \right)}} - {{{z}}^{(k + 1)}}} \right) \end{aligned} \right\}$$ (18)

      其中,上标$k$表示第$k$次迭代,式(18)中前两个子问题可以通过对函数${L_\rho }({{x}},{{z}},{{\alpha}} )$中的${{x}}$${{z}}$分别求导并令导数为0求得,经求解式(18)中${{x}},{{z}}$以及${{\alpha}} $的迭代式为

      $$ \left.\begin{split}{{x}}^{\left(k+1\right)}=&\bigg[{\left({{F}}^{\left(3\right)}\otimes {{F}}^{(2)}\otimes {{F}}^{\left(1\right)}\right)}^{\rm{H}}\\ &\cdot{\left({{F}}^{\left(3\right)}\otimes {{F}}^{(2)}\otimes {{F}}^{\left(1\right)}\right)+\rho {{I}}_{UVW}\bigg]}^{-1}\\ &\cdot \left[{\left({{F}}^{\left(3\right)}\otimes {{F}}^{(2)}\otimes {{F}}^{\left(1\right)}\right)}^{\rm{H}}{{s}}+\rho {{z}}^{(k)}-{{{\alpha}} }^{(k)}\right]\\ {{z}}^{(k+1)}=&{\rm{ST}}\left({{x}}^{\left(k+1\right)}+\frac{{{{\alpha}} }^{(k)}}{\rho };\frac{\lambda }{\rho }\right)\\ {{{\alpha}} }^{(k+1)}=&{{{\alpha}} }^{(k)}+\rho \left({{x}}^{\left(k+1\right)}-{{z}}^{(k+1)}\right) \end{split}\!\!\right\}$$ (19)

      其中,ST为软阈值(Soft Threshold)函数,其表达式为${\rm{ST}}\left( {x,a} \right) = \left( {x/\left| x \right|} \right)\max\left( {\left| x \right| - x,0} \right)$。将F(1)=T1·F1, ${{{F}}^{(2)}} = {{{T}}_2}{{{F}}_2}$, ${{{F}}^{(3)}} = {{{T}}_3}{{{F}}_3}$代入式(19)可得

      $$ \begin{split} {{{x}}^{\left( {k + 1} \right)}} = &\left[ {{{\left( {{{{F}}_3} \otimes {{{F}}_2} \otimes {{{F}}_1}} \right)}^{\rm{H}}}\left( {{{{B}}_3} \otimes {{{B}}_2} \otimes {{{B}}_1}} \right)} \right. \\ &\cdot{ {\left( {{{{F}}_3} \!\otimes\! {{{F}}_2} \!\otimes\! {{F}}} \right) \!+ \!\rho {{{I}}_{UVW}}} \Big]^{{\rm{ - }}1}} \! \cdot \! \! \left[ {{{\left( {{{{F}}_3}\! \otimes\! {{{F}}_2} \otimes {{F}}} \right)}^{\rm{H}}}} \right. \\ &\cdot\left. {{{\left( {{{{T}}_3} \otimes {{{T}}_2} \otimes {{{T}}_1}} \right)}^{\rm{H}}}{{s}} + \rho {{{z}}^{(k)}} - {{{\alpha}} ^{(k)}}{\kern 1pt} } \right]\\[-15pt]\end{split} $$ (20)

      其中,${{{B}}_1} = {{T}}_1^{\rm{H}}{{{T}}_1},{{{B}}_2} = {{T}}_2^{\rm{H}}{{{T}}_2},{{{B}}_3} = {{T}}_3^{\rm{H}}{{{T}}_3}$,通过化简可得

      $$\begin{split} {{{x}}^{\left( {k + 1} \right)}} = &{\left( {{{{F}}_3} \otimes {{{F}}_2} \otimes {{{F}}_1}} \right)^{\rm{H}}}\left[ {\left( {{{{B}}_3} \otimes {{{B}}_2} \otimes {{{B}}_1}} \right)} \right. \\ &+{ { \rho {{{I}}_{UVW}}} \Big]^{{\rm{ - }}1}}{\kern 1pt} \left[ {{{\left( {{{{T}}_3} \otimes {{{T}}_2} \otimes {{{T}}_1}} \right)}^{\rm{H}}}{{s}}} \right. \\ &+ \left( {{{{F}}_3} \otimes {{{F}}_2} \otimes {{{F}}_1}} \right) \left. {\left( {\rho {{{z}}^{(k)}} - {{{\alpha}} ^{(k)}}} \right)} \right] \end{split} $$ (21)

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

      $$\begin{split} {{\cal{X}}^{\left( {k + 1} \right)}} = &\Big\{ {\Big[ {\left( {{\cal{S}}{ \times _1}{{{T}}_1}{ \times _2}{{{T}}_2}{ \times _3}{{{T}}_3}} \right)} } \\ &\left. { + \; \left( {\left( {\rho {{\cal{Z}}^{(k)}} - {{\cal{A}}^{(k)}}} \right){ \times _1}{{{F}}_1}{ \times _2}{{{F}}_2}{ \times _3}{{{F}}_3}} \right)} \right] \\ & {{ \times _1}{{F}}_1^{\rm{H}}{ \times _2}{{F}}_2^{\rm{H}}{ \times _3}{{F}}_3^{\rm{H}}} \Big\} \oslash \left[ {{\cal{G}} + \rho {{\bf{1}}_{U \times V \times W}}} \right] \end{split} $$ (22)

      其中,${{\bf{1}}_{U \times V \times W}}$表示元素全为1、维数为$U \times V \times W$的三维张量,$ \oslash $表示张量的元素除法,${\cal{G}}$表示对接收回波三维方向的采样,其值设为0或1,分别表示是否被采样到。同理式(19)中${{{z}}^{(k + 1)}}$${{{\alpha}} ^{(k + 1)}}$的迭代式同样可以写成如式(23)和式(24)的张量形式

      $$ {{\cal{Z}}}^{(k+1)}={\rm{ST}}\left({{\cal{X}}}^{\left(k+1\right)}+\frac{{{{\alpha}} }^{(k)}}{\rho };\frac{\lambda }{\rho }\right)$$ (23)
      $${{\cal{A}}^{(k + 1)}} = {{\cal{A}}^{(k)}} + \rho \left( {{{\cal{X}}^{\left( {k + 1} \right)}} - {{\cal{Z}}^{(k + 1)}}} \right)$$ (24)

      联合迭代式(22)—式(24),就可得到图像${\cal{X}}$。初始参数设置如下:${\cal{Z}}$${\cal{A}}$的初值设定为${\bf{0}}$$\lambda $的值根据数据进行调整,$\rho $的值取1。

      将三维形式进一步推广,假设多维张量的维数为${U_1} \times {U_2} \times \cdots \times {U_i}$,则张量${{\cal{X}}^{\left( {k + 1} \right)}}$的更新表达式为

      $$ \left. \begin{aligned} &{{\cal{X}}}^{\left(k+1\right)}=\bigg\{\bigg[\left({\cal{S}}{\times }_{1}{{T}}_{1}{\times }_{2}{{T}}_{2}\times\cdots {\times }_{i}{{T}}_{i}\right)\\ &\quad +\left(\left(\rho {{\cal{Z}}}^{(k)}-{{\cal{A}}}^{(k)}\right){\times }_{1}{{F}}_{1}{\times }_{2}{{F}}_{2}\times\cdots {\times }_{i}{{F}}_{i}\right)\bigg]{\times }_{1}\\ &{{F}}_{1}^{\rm{H}}{\times }_{2}{{F}}_{2}^{\rm{H}}\times\cdots {\times }_{i}{{F}}_{i}^{\rm{H}}\bigg\}\oslash \left({\cal{G}}+\rho {\bf{1}}_{{U}_{1}\times {U}_{2}\times\cdots \times {U}_{i}}\right)\\ &{{\cal{Z}}}^{(k+1)}={\rm{ST}}\left({{\cal{X}}}^{\left(k+1\right)}+\frac{{{{\alpha}} }^{(k)}}{\rho };\frac{\lambda }{\rho }\right)\\ &{{\cal{A}}}^{(k+1)}={{\cal{A}}}^{(k)}+\rho \left({{\cal{X}}}^{\left(k+1\right)}-{{\cal{Z}}}^{(k+1)}\right) \end{aligned} \right\}$$ (25)
    • MD-ADMM算法中图像${\cal{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},其中L2L3分别为第1层循环和第2层循环的迭代次数。根据文献[5],算法MD-SL0的计算复杂度为O[(MNKU+UVNK+UVKW+UVWM+MNVW+MNKW)L2L3]。在本实验中,迭代次数${L_2}{L_3} \gg {L_1}$,且算法DR-2D-SL0的单次迭代计算复杂度最高,因此3种算法的计算复杂度排序为:MD-ADMM<MD-SL0<DR-2D-SL0。

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

      图  4  仿真目标三维散点图

      Figure 4.  3D scatter of simulation target

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

      图  5  完整回波数据图像三视图

      Figure 5.  Three views of image with the complete echo

      图  6  回波采样形式

      Figure 6.  Undersampling masks of random sampling and block sampling

      首先采用随机采样方式,对稀疏度(每个维度的稀疏度相同)分别为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都得到了质量较高的图像。

      图  7  稀疏度为50%时图像

      Figure 7.  Image when sparsity is 50%

      图  8  稀疏度为33.3%时的图像

      Figure 8.  Image when sparsity is 33.3%

      图  9  稀疏度为25%时的图像

      Figure 9.  Image when sparsity is 25%

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

      表 1  随机稀疏采样条件下数值结果

      Table 1.  Numerical results under random sparse sampling condition

      稀疏度算法图像熵PSNR计算时间
      50.0%RD9.49034.0750.012
      DR-2D-SL03.13857.71435.053
      MD-SL03.13857.7146.954
      MD-ADMM3.00661.6783.518
      33.3%RD10.50528.4300.010
      DR-2D-SL03.20052.70229.091
      MD-SL03.20052.70212.546
      MD-ADMM2.95953.8886.088
      25.0%RD10.96526.0030.0131
      DR-2D-SL03.37348.75554.008
      MD-SL03.37348.75528.509
      MD-ADMM3.02349.38611.996
      $$E({\cal{X}}){\rm{ = }} - {\rm{Sum}}\left\{ {\frac{{{{\cal{X}}^2}}}{{\rm{C}}} \odot \ln\left[ {\frac{{{{\cal{X}}^2}}}{{\rm{C}}}} \right]} \right\}$$ (26)

      其中,${\rm{Sum}}( \cdot )$表示张量的所有元素之和,$C = \displaystyle\sum\nolimits_{u = 1}^U {\displaystyle\sum\nolimits_{v = 1}^V {\displaystyle\sum\nolimits_{w = 1}^W {x_{uvw}^2} } }$, $U,V,W$分别为张量${\cal{X}}$的维数。假设对${\cal{X}}$进行了归一化,则${\cal{X}}$相对于参考图像${\cal{H}}$的均方误差可以表示为

      $${\rm{MSE}} = \frac{1}{{UVW}}\sum\limits_{u = 1}^U {\sum\limits_{v = 1}^V {\sum\limits_{w = 1}^W {\left\| {{x_{uvw}} - {h_{uvw}}} \right\|} } } $$ (27)

      其中,${h_{uvw}}$为参考三维图像${\cal{H}}$的元素。定义PSNR为

      $${\rm{PSNR}} = 10 \cdot {\lg}\left( {{{{1^3}} / {{\rm{MSE}}}}} \right)$$ (28)

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

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

      表 2  块稀疏采样条件下数值结果

      Table 2.  Numerical results under block sparse sampling condition

      稀疏度算法图像熵PSNR计算时间
      50.0%RD8.33632.8850.012
      DR-2D-SL03.24651.23241.907
      MD-SL03.24651.2328.029
      MD-ADMM3.03753.2743.660
      33.3%RD8.55327.4950.0098
      DR-2D-SL03.70042.18842.518
      MD-SL03.70042.18817.954
      MD-ADMM3.23342.3365.758
      25.0%RD9.50026.9060.0092
      DR-2D-SL03.89445.19651.463
      MD-SL03.89445.19625.702
      MD-ADMM3.33747.20111.536

      图  10  稀疏度为50%的块稀疏采样图像

      Figure 10.  The image of a sparsity of 50% by block sampling

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

      表 3  不同信噪比条件下数值结果

      Table 3.  Numerical results under different SNR conditions

      信噪比算法图像熵PSNR计算时间
      –5 dBRD11.804218.3200.0177
      DR-2D-SL07.321837.63949.998
      MD-SL07.321837.63923.540
      MD-ADMM4.77943.35210.775
      0 dBRD11.53523.4520.013
      DR-2D-SL05.408643.28548.776
      MD-SL05.408643.28523.303
      MD-ADMM2.97747.44910.715
      10 dBRD11.11627.1450.0124
      DR-2D-SL03.94845.61948.754
      MD-SL03.94845.61923.484
      MD-ADMM3.06647.56110.521

      图  11  稀疏度为25%信噪比为–5 dB的图像

      Figure 11.  The image when sparsity is 25% and SNR=–5 dB

      图  12  稀疏度为25%信噪比为0 dB图像

      Figure 12.  The image when sparsity is 25% and SNR=0 dB

      图  13  稀疏度为25%信噪比为10 dB的图像

      Figure 13.  The image when sparsity is 25% and SNR=10 dB

    • MIMO雷达实验系统目前还在搭建中,还存在发射信号带宽过窄、收发阵列同步性差等问题,导致MIMO-ISAR回波目前还难以获取,也是下一步要着重解决的问题。因此采用二维Yak-42飞机ISAR实测数据,以验证所提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算法基本失效,算法MD-SL0, MD-ADMM在二维稀疏采样条件下依然能够得到清晰图像,但相比MD-SL0算法,MD-ADMM算法得到的图像熵最小,峰值信噪比最大,并且计算时间最短,进一步验证了所提算法的有效性。

      表 4  实测数据不同信噪比条件下的数值结果

      Table 4.  Numerical results of measured data for different signal-to-noise ratio conditions

      信噪比算法图像熵PSNR计算时间
      RD9.63929.8160.010
      MD-SL06.45940.37810.553
      MD-ADMM5.29641.3643.743
      0 dBRD10.24226.9090.0131
      MD-SL07.57135.24112.064
      MD-ADMM6.18339.0904.274

      图  14  实测数据结果

      Figure 14.  Measured ISAR data results

      图  15  信噪比为0 dB实测数据结果

      Figure 15.  Measured ISAR data results under SNR=0 dB

    • 本文提出一种多维ADMM稀疏恢复算法,本算法可以用于恢复多维稀疏信号,而无需将多维信号转换为一维,这极大地减少了存储量和计算量。通过该算法,实现了一种实用的MIMO-ISAR成像。与其他算法相比,本算法具有存储容量小、计算效率高、成像质量好的优点。仿真和实测数据结果均证明了本算法的有效性。

参考文献 (15)

目录

    /

    返回文章
    返回