1.   引言

在合成孔径雷达(Synthetic Aperture Radar, SAR)研究领域，针对特定目标的自动分类和识别一直以来都是最为热门的研究方向之一。如何稳健并高效地获取目标相关有用特征，对实现目标分类和识别至关重要。近年来，由Donoho等人^[1,2]提出的压缩感知(Compressed Sensing, CS)技术，为SAR成像稀疏特征增强提供了有效的技术实现途径^[3-8]。压缩感知技术可有效降低高分辨成像对系统数据采样率的要求，从而简化系统复杂度，同时降低在实际雷达系统应用中，由于功能切换、系统误差以及电子干扰等因素造成的影响，为后续目标分类识别提供便利条件。

2007年，Baraniuk^[3]首次成功应用CS理论实现高分辨SAR稀疏成像，时至今日，CS已应用于SAR^[4,5]、SAR地面运动目标成像(Ground Moving Target Imaging, GMTIm)^[6,7]以及逆SAR(Inverse SAR, ISAR)^[8]多种模式。近十年发展中，相关稀疏特征增强成像算法可大致分为以下3类：

(1) 贪婪(Greedy)类算法：该类算法主要通过迭代搜索均方根误差最小的局部稀疏最优解，并在解的稀疏度小于某一经验值时停止迭代。正交匹配追踪(Orthogonal Matching Pursuit, OMP)算法是最为经典的贪婪类算法之一^[9]，文献[4]将贪婪算法应用在宽测绘带的雷达信号处理中, 验证了此方法在成像效率上的优势。尽管该类算法计算复杂度低，运算速度快，但解的稀疏度并不高，因此针对稀疏SAR特征增强的精度有限。

(2) 凸优化类算法：该类算法通过最小化目标的${\ell _1}$范数，结合基追踪噪声抑制(Basis Pursuit De-Noising, BPDN)，实现对稀疏信号的有效恢复。数学推导已经证明^[10]，${\ell _1}$范数可认为是最优的稀疏近似解。凸优化类算法由于应用了${\ell _1}$范数正则化模型，在稀疏信号恢复精度上明显优于贪婪算法，并在SAR领域广泛应用。文献[5]对SAR成像中的相位误差问题使用该类算法，实现了自聚焦。然而，经典的凸优化(ConVeX, CVX)主要利用了递归最小二乘的思想，在求解${\ell _1}$范数正则化模型时每步迭代都要进行高维矩阵的求逆运算，因此存在计算效率低，对系统误差敏感等问题。

(3) 贝叶斯类算法：该类算法基于非监督的统计机器学习理论，通过引入SAR图像目标稀疏先验分布，利用分层模型进行贝叶斯推理，从而重建稀疏特征的后验概率密度函数估计，最终获得高精度的稀疏信号恢复性能。文献[6,7]将贝叶斯压缩感知(Bayesian Compress Sensing, BCS)方法应用于SAR-GMTIm模式，得到了稀疏特征增强后的地面运动目标图像。文献[8]通过贝叶斯分层模型，实现了ISAR成像稀疏特征增强，同时对平动误差进行了统计建模从而重建目标的聚焦特征。BCS稀疏特征增强方法虽然精度较高，但是在求解过程中难以避免求逆运算，因此存在着计算复杂度大，运算效率偏低等问题。

综上所述，有必要研究一种能在SAR稀疏特征增强中，获得精度和效率之间良好平衡的方法。2011年，Boyd等人^[11]说明了交替方向多乘子方法(Alternating Direction Method of Multipliers, ADMM)可适用于大规模分布式计算系统及优化问题。这种方法最早是由Gabay和Mercier于1976年提出的，由于当时相关大规模计算系统的技术有限，所以ADMM在相当长的时间内并未得到重视。近几年随着高维大数据的出现与相关计算技术与条件的成熟，ADMM在计算效率上的优势逐渐体现出来。ADMM通过对优化变量进行对偶分解，并结合增广的拉格朗日多乘子方法建立优化方程，利用高斯-赛德尔(Gaussian-Seidel)思想进行交替方向迭代，将复杂的优化问题分解为多个简单的子优化问题，通过分解-调和(Decomposition-Coordination)方式逼近最优解。文献[11]中严格分析了ADMM的迭代残差及算法收敛性，可以保证算法在迭代运算中严格收敛，更加稳健快速地达到收敛条件，因此ADMM在SAR稀疏信号增强方面具有优异的稳健性。此方法目前已被广泛应用在图像处理领域，但是相关实验大多都局限在实数图像域，由于SAR成像回波数据的复数特性，常规的ADMM方法无法直接对复数SAR回波数据的稀疏成像。

为了便于检验算法在目标稀疏特征增强中的性能，本文推导得出了SAR, SAR-GMTIm以及ISAR成像模式的通用回波信号解析式，并基于此应用复数ADMM(Complex ADMM, C-ADMM)的框架进行求解。C-ADMM首先对包含保真项和增广拉格朗日项的岭回归(ridge regression)问题求解，然后由${\ell _1}$范数正则项和拉格朗日项推导复数${\ell _1}$范数邻近算子，即复数软阈值^[12]，最后更新岭回归解和软阈值联合最小化的对偶变量，从而实现对复数SAR数据的稀疏特征增强。本文实验应用多种SAR模式实测数据，通过对比C-ADMM与CVX, BCS算法的稀疏成像结果，定量和定性地分析了C-ADMM在性能上的优势，验证了本文所提算法的稳健性、高效性和信号模型的通用性。

2.   回波信号模型

本文提出SAR, SAR-GMTIm以及ISAR模式通用的回波信号模型，可实现多种模式回波数据的稀疏特征增强。如图1所示，图中分别给出了SAR, SAR-GMTIm以及ISAR成像的几何模型。根据合成孔径原理，3种模式的成像机理都是通过载台或者目标的运动形成虚拟合成孔径并最终实现高分辨成像。其中，SAR成像通过载机平台沿预定航线以速度$v$飞行，地面静止目标场景中心${O'}$与载机航线参考距离矢量为${{{R}}_0}$。为了进一步适应地面运动目标的成像需求，特采用多通道天线，以在成像处理过程中消除地面静止场景回波，即杂波的影响，从而便于检测运动目标^[13]。设${{{q}}_i}$为原点到第$i$通道天线相位中心(Antenna Phase Center, APC)位置矢量，${{{q}}_0}$为参考通道位置矢量，如图1所示。假设地面目标位置矢量为${{{r}}_{\rm t}}\left( t \right)$，待成像目标与天线参考通道APC在$t$时刻的距离可以表示为

图 1  SAR成像几何模型

下载: 全尺寸图片 幻灯片

$$R\left( t \right) = \left| {{{r}_{\rm t}}\left( t \right) - {{q}_0}\left( t \right)} \right|$ $

(1)

其中，$t$表示方位时间变量。

(1) SAR模式

在SAR模式下，地面成像场景静止，此时实时目标斜距可按照泰勒公式在参考斜距$\left| {{{R}_0} - {{q}_0}\left( t \right)} \right|$处展开

$$\begin{split} R\left( t \right) =& \left| {{{R}_0} + {{r}_p} - {{q}_0}\left( t \right)} \right|\\ =& \left| {{{R}_0} - {{q}_0}\left( t \right)} \right| + {{r}_p}\frac{{{{R}_0} - {{q}_0}\left( t \right)}}{{\left| {{{R}_0} - {{q}_0}\left( t \right)} \right|}} + o\left( t \right) \end{split}$ $

(2)

其中，${{r}_p}$=$\left( {{x_p},{r_p}} \right)$表示第$p$个地面静止散射点在方位向和距离向的位置矢量，$o\left( t \right)$表示在参考斜距$\left| {{{R}_0} - {{q}_0}\left( t \right)} \right|$处的2阶及以上展开项，在高分辨小场景成像条件下高阶项的影响并不大，可以忽略。根据公式${{R}_0}$=${R_{\rm{0}}}{r}$, ${{q}_0}\left( t \right)$=$u{x}$，其中$u$=$vt$, ${r}$和${x}$分别表示距离向和方位向单位矢量，可以得到静止场景的SAR多点累积观测回波

$$\begin{split}{S_0}\left( {k,t} \right) =& P\left( k \right)\sum\limits_{p = 1}^P \exp \Biggr\{ - {\rm{j}}k\Biggr( \sqrt {{R_0}^2 + {u^2}} \\ & + \frac{{{R_0}{r_p} - u{x_p}}}{{\sqrt {{R_0}^2 + {u^2}} }}\Biggr) \Biggr\} \end{split} $ $

(3)

式(3)经过方位向dechirp操作及极坐标插值算法(Polar Formation Algorithm, PFA)^[7]处理，令${k_r}$=${{{R_0}k}/{\sqrt {{R_0}^2 + {u^2}} }}$, ${k_x}$=${{ - uk}/{\sqrt {{R_0}^2 + {u^2}} }}$，可以得到回波2维波数域表达式为

$${S_0}({k_r},{k_x}) = \sum\limits_{p = 1}^P {\exp \left( { - {\rm{j}}{k_x}{x_p} - {\rm{j}}{k_r}{r_p}} \right)} $ $

(4)

在PFA处理中有${k_x}$=${{ - v{t'}{k_0}}/{{R_0}}}$，其中$t'$表示插值后的时域变量，与${k_x}$线性相关。此时，可得SAR模式的数据域表达式为

$$\begin{split} {S_0}(\hat r,t') &= \sum\limits_{p = 1}^P {{\rm{sinc}}\left( {\hat r - {r_p}} \right)\exp \left( { - {\rm{j}}{k_x}{x_p}} \right)} \\ &= \sum\limits_{p = 1}^P {{\rm{sinc}}\left( {\hat r - {r_p}} \right)\exp \left( {{\rm{j2}}\pi \frac{{2v{x_p}}}{{\lambda {R_0}}}t'} \right)} \end{split}$ $

(5)

其中，回波信号可以表示为$P$个地面散射点回波信号的叠加求和形式，其中${\rm {sinc}}$函数为第$p$个散射点的距离向包络，$\exp$函数为方位线性相位项。

(2) SAR-GMTIm模式

在SAR-GMTIm模式下，需要考虑地面动目标与载机平台的相对运动，平台与地面第$k$个运动目标的斜距历程可以写为

$$\begin{split} R\left( t \right)\!=\! & \left| {{{R}_0} \!+ \!{{r}_k} \!+\! {{v}_k}t \!+\! \frac{1}{2}{{a}_k}{t^2} \!-\! {{q}_0}\left( t \right)} \right|\\ =& \left| {{{R}_0} \!-\! {{q}_0}\left( t \right)} \right| \!+\! \left( {{{r}_k} \!+\! {v}_k^{}t \!+\! \frac{1}{2}{{a}_k}{t^2}} \right)\frac{{{{R}_0} \!-\! {{q}_0}\left( t \right)}}{{\left| {{{R}_0} \!-\! {{q}_0}\left( t \right)} \right|}} \\ &+ \frac{{{{\left( {{{r}_k} + {{v}_k}t + \frac{1}{2}{{a}_k}{t^2}} \right)}^2}}}{{2\left| {{{R}_0} - {{q}_0}\left( t \right)} \right|}} + o\left( t \right)\\[-15pt] \end{split}$ $

(6)

其中，${{r}_k}$,${{v}_k}$和${{a}_k}$分别为目标的初始位置、速度和加速度矢量，可分解为

$$\left. \begin{array}{l} {{r}_k} = x_k^{}{x} + r_k^{}{r}\\ {{v}_k} = v_k^x{x} + v_k^r{r}\\ {{a}_k} = a_k^x{x} + a_k^r{r} \end{array} \right\}$ $

(7)

由于本文GMTIm模式采用沿航向多通道SAR体制，根据通道间地面静止目标回波(杂波)的方位不变性，可实现杂波抑制，从而有效提高对运动目标成像的信噪比。式(7)经过对每个通道的SAR回波进行成像处理，再利用DPCA方法^[14]，将任意通道与参考通道进行相消，便可以得到杂波抑制后总回波为

$$ \begin{split} {S_{i {\scriptsize{-}} 0}}\left( {\hat r,t'} \right) =& \sum\limits_{k = 1}^K {{A_k}} {\rm{sinc}}\left( {\hat r - {r_k}} \right)\exp \left( {{\rm{j}}\frac{{\psi _i^{}}}{{\rm{2}}}} \right)\\ &\cdot\exp \left( {{\rm{j}}2\pi f_d^kt' + {\rm{j}}\pi \gamma _d^k{t^{'2}}} \right) + {\rm{C}}{{\rm{N}}_{i - 0}}\left( {\hat r,t'} \right) \end{split}$ $

(8)

其中，${S_{i {\tiny{-}} 0}}\left( {\hat r,{t'}} \right)$可以分解为由地面静止目标产生的残留杂波${\rm{C}}{{\rm{N}}_{i - 0}}\left( {\hat r,{t'}} \right)$与其余由$K$个观测目标产生的运动回波之和，第$k$个运动目标的回波可以表示为距离向包络${\rm{sinc}}$函数与运动引起的方位向的线性相位乘积的形式，其中${\psi _i}$=$\left( {{\rm{4}}\pi {\rm{/}}{\lambda _0}} \right)\left( {{v_r}/v} \right){d_i}$表示干涉相位，${d_i}$为${{q}_{\rm{i}}}$与${{q}_{\rm{0}}}$通道之间的距离。${\lambda _0}$表示工作波长，其中$f_d^k$和$\gamma _d^k$分别为目标运动造成的方位多普勒频率和调频率，即

$$\left. {\begin{array}{*{20}{l}} {f_d^k \!=\! \dfrac{{2v}}{{{\lambda _0}}}\dfrac{{{x_k}}}{{{R_0}}} \!-\! \dfrac{2}{{{\lambda _0}}}\dfrac{{\left( {{R_0} + {r_k}} \right)v_k^r}}{{{R_0}}}}\\ {\gamma _d^k \!=\! \dfrac{2}{{{\lambda _0}}}\dfrac{{{v^2}}}{{{R_0}}} \!-\! \dfrac{2}{{{\lambda _0}}}\dfrac{{{{\left( {v - v_k^x} \right)}^2}}}{{{R_0}}} \!-\! \dfrac{2}{{{\lambda _0}}}\dfrac{{\left( {{R_0} \!+\! {r_k}} \right)a_k^r}}{{{R_0}}}} \end{array}} \right\}$ $

(9)

其中，多普勒频率和调频率可通过吕氏分布($\rm LV's$ Distribution, LVD)时频表示方法进行估计^[15]。

(3) ISAR模式

合成孔径概念可以推广至空中运动目标成像中，ISAR成像几何模型如图1(b)所示，经过补偿平动误差分量后^[16]，目标可等效为转台模型，并以均匀角速度$\omega$转动。在小转角假设下，目标在相干积累时间内的转角$\theta \left( t \right) \approx 0$，可得目标散射点与雷达距离历程为

$$\begin{split}R\left( t \right) =& {R_0} + {r_p}\cos \theta \left( t \right) + {x_p}\sin \theta \left( t \right) \\ \approx & {R_0} + {r_p} + {x_p}\omega t\end{split}$ $

(10)

此时与SAR模式的处理流程相同，可得到回波

$$\begin{split} S(\hat r,t) =& \sum\limits_{p = 1}^P {{\rm{sinc}}\left( {\hat r - {r_p}} \right)\exp \left( { - {\rm{j}}k{x_p}\omega t} \right)} \\ \approx &\sum\limits_{p = 1}^P {{\rm{sinc}}\left( {\hat r - {r_p}} \right)\exp \left( {{\rm{j2\pi }}\frac{{2{x_p}\omega }}{\lambda }t} \right)}\end{split} $ $

(11)

由式(11)可见，ISAR与SAR模式类似，相应回波也可表示为距离向${\rm{sinc}}$包络与方位向线性相位乘积的$P$点累积求和的形式。

综上所述，多种SAR成像模式的距离压缩域数据$S\left( {\hat r,t} \right)$或$S\left( {\hat r,{t'}} \right)$，均可表示为距离向${\rm{sinc}}$包络与方位向线性相位乘积，并进行多点累积求和的形式，因此SAR回波成像类似于线性回归模型

$${Y} = {AX} + {{C}}_{\rm{N}}$ $

(12)

其中，${Y}$表示距离压缩预处理后的SAR数据，即推导得出的$S\left( {\hat r,t} \right)$或$S\left( {\hat r,{t'}} \right)$, ${X}$为待恢复的目标或图像，相应像素对应了方位向和距离向位置，${{C}_{\rm{N}}}$为加性的噪声、杂波和干扰等，${A}$为方位向傅里叶字典，根据上述讨论字典${A}$可表示为

$$\left. {\begin{array}{*{20}{l}} {{A} \!=\! {A}\left( {{{\gamma }_d}} \right) = \left[ {{{A}_1},{{A}_2}, ··· ,{{A}_k}} \right]}\\ {{{A}_k} \!=\! {{A}_0} \odot {B}\left( {{\gamma _d}^k} \right)}\\ {{{A}_0} \!=\! \left[ {{a}\left( {{f_d}\left( 1 \right)} \right),{a}\left( {{f_d}\left( 2 \right)} \right), ··· ,{a}\left( {{f_d}\left( n \right)} \right)} \right]}\\ {a}\left( {{f_d}\left( n \right)} \right) = \left[ {{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}{f_d}\left( n \right){t_1}}},\right.\\ \left.{{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}{f_d}\left( n \right){t_2}}}, ··· ,{{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}{f_d}\left( n \right){t_n}}} \right]^{\rm{T}}\\ {{B}\left( {{\gamma _d}^k} \right) \!=\! {{\left[ {{b}\left( {{\gamma _d}^k} \right),{b}\left( {{\gamma _d}^k} \right), ··· ,{b}\left( {{\gamma _d}^k} \right)} \right]}_{1 \times N}}}\\ {{b}\left( {{\gamma _d}^k} \right) \!=\! {{\left[ {{{\rm{e}}^{{\rm{j\pi }}{\gamma _d}^k{t_1}^2}},{{\rm{e}}^{{\rm{j\pi }}{\gamma _d}^k{t_2}^2}}, ··· ,{{\rm{e}}^{{\rm{j\pi }}{\gamma _d}^k{t_N}^2}}} \right]}^{\rm{T}}}} \end{array}}\!\!\! \right\}\!\!\!\!\!\!$ $

(13)

其中，${\left[ \cdot \right]^{\rm{T}}}$表示矩阵转置操作符，$\odot$为Hadamard积操作符。在SAR和ISAR成像模式下，${A}$为常规的傅里叶字典，即${A}$=${{A}_0}$，在GMTIm的成像模式下，${A}$为2阶参数化傅里叶字典，即${A}$=${A}\left( {{\gamma _d}} \right)$ 。

3.   C-ADMM

由于观测回波噪声的影响，直接求解如式(12)所示的SAR成像通用回波模型，如${X} = {{A}^{ - 1}}{Y}$，将出现严重的误差，${\left[ \cdot \right]^{{\rm{ - 1}}}}$ 表示矩阵求逆运算。因此Tibshirani^[17]在1996年提出了最小绝对收缩和变量选择算子(the Least Absolute Shrinkage and Selection Operator, LASSO)模型，在原来只有保真项的基础上增加了${\ell _1}$范数正则项，更有利于选择理想的解，LASSO模型为

$$\hat{ X} = \arg \mathop {\min }\limits_{{X}} \left( {\frac{1}{2}\left\| {{Y} - {AX}} \right\|_{\rm{F}}^2{\rm{ + }}\lambda {{\left\| {X} \right\|}_1}} \right)$ $

(14)

其中，${\left\| \cdot \right\|_1}$和${\left\| \cdot \right\|_{\rm{F}}}$分别表示${\ell _1}$范数和Fibonacci范数，$\lambda$为正则项系数，通常$\lambda > 0$。式(14)中第1项为保真项，表示恢复后雷达图像和原始图像的逼近程度，第2项为正则项，表征了SAR图像中强散射点具有的稀疏先验信息。通过调整参数$\lambda$可调节稀疏解的稀疏度，同时达到较好的抑噪性能。Boyd等人^[11]在2012年较为系统地提出了基于ADMM方法求解LASSO问题的操作步骤，同时结合了对偶分解和增广的拉格朗日多乘子方法，将优化问题分步进行求解，从而在提高计算效率的同时获得稳健的收敛性能。ADMM方法的优势在于将一个复杂的优化问题分解为若干个子问题，在实现分布式优化的同时可以调和全局性能，从而达到对最优解的逼近。依据ADMM算法理论，LASSO可以分解为两个子问题$f\left( {X} \right)$=$\left( {1/2} \right)\left\| {{Y} - {AX}} \right\|_{\rm{F}}^2$和$g\left( {Z} \right)$=$\lambda {\left\| {Z} \right\|_1}$，此时，式(14)的优化问题可以等效表示为

$${\rm{minimize}}\;f\left( {X} \right) + g\left( {Z} \right),\quad {\rm{s}}.{\rm{t}}.\;{X} - {Z} = 0$ $

(15)

由此，建立带有增广拉格朗日项的优化方程，如

$$ \begin{split} {{{L}}_\rho }\left( {{X},{Z},{U}} \right){\rm{ = }}&f\left( {X} \right){\rm{ + }}g\left( {Z} \right){\rm{ + }}{{U}^{\rm{H}}}\left( {{X} - {Z}} \right)\\ &{\rm{ + }}\left( {\rho /2} \right)\left\| {{X} - {Z}} \right\|_{\rm{F}}^2\end{split}$ $

(16)

其中，${\left[ \cdot \right]^{\rm{H}}}$表示矩阵共轭转置操作符，$\rho$为拉格朗日乘子系数，引入的拉格朗日乘子项为优化方程式(17)加上了等式约束，将目标解限定在相应的可行域之内，提升了求解效率，即算法的稳健性。此时，求得ADMM的解可表示为

$$ \!\!\!\left. {\begin{array}{*{20}{l}} {{{X}^{k \!+\! 1}} \!= \!{\rm{arg}}\mathop {{\rm{min}}}\limits_{{X}} \left\{ {f\left( {X} \right) \!+\! \dfrac{\rho }{2}\left\| {{X} \!-\! {{Z}^k} \!+\! {{U}^k}} \right\|_{\rm{F}}^2} \right\}}\\ {{{Z}^{k \!+\! 1}} \!=\! {\rm{arg}}\mathop {{\rm{min}}}\limits_{{Z}} \left\{ {g\left( {Z} \right) \!+\! \dfrac{\rho }{2}\left\| {{{X}^{k \!+\! 1}} \!-\! {Z} \!+\! {{U}^k}} \right\|_{\rm{F}}^2} \right\}}\\ {{{U}^{k \!+\! 1}} \!=\! {{U}^k} \!+\! {{X}^{k + 1}} \!-\! {{Z}^{k \!+\! 1}}} \end{array}} \!\!\!\!\right\}\!\!$ $

(17)

其中，上标$k$表示矩阵第$k$次迭代值。变量${X}$和${Z}$通过交替方向的迭代更新可达到联合最小化，而变量${U}$则利用了${X}$和${Z}$的迭代值获得对偶变量更新，在求解过程中，拉格朗日乘子分别作用于变量${X}$和${Z}$，然后通过${X}$和${Z}$的迭代间接作用于变量${U}$，从而提高了收敛速率与稳健性。下面分别针对式(17)的3个步骤说明如下：

步骤 1　岭回归 由式(17)中第1个等式，变量${X}$的优化问题等同于保真项与拉格朗日项两个Fibonacci范数的联合最小化问题，该最小化问题可理解为岭回归问题，即

$$\begin{split}{{X}^{k + 1}}=&{\rm{ arg}}\mathop {{\rm{min}}}\limits_{{X}} \left\{ \frac{1}{2}{\rm{Tr}}\left[ {{{\left( {{Y} - {AX}} \right)}^{\rm{H}}}\left( {{Y} - {AX}} \right)} \right] \right.\\ & + \frac{\rho }{2}{\rm{Tr}}\left[ {{{\left( {{X} - {{Z}^k} + {{U}^k}} \right)}^{\rm{H}}}\left( {{X} - {{Z}^k} + {{U}^k}} \right)} \right] \} \end{split}$ $

(18)

其中${\rm{Tr}}\left[ \cdot \right]$表示求迹运算符。由式(18)可见，最小化目标函数为凸函数，因此最优值可通过对目标函数求导得到，解得

$${{X}^{k + 1}} \!=\! {\left( {{{A}^{\rm{H}}}{A} \!+\! \rho {I}} \right)^{ - 1}}\left\{ {{{A}^{\rm{H}}}{Y} \!+\! \rho \left( {{{Z}^k} \!-\! {{U}^k}} \right)} \right\}$ $

(19)

步骤 2　${\ell _1}$范数正则化 再由式(17)，变量${Z}$的优化问题可表示求解${\ell _1}$范数正则化问题，区别于式(14)的LASSO问题，这里的正则化面向的是增广的拉格朗日项。由${\ell _1}$范数正则化，可实现对成像场景稀疏特征的先验表征，利用式(18)求解${X}$中的梯度下降思想，得

$$\begin{split} {{Z}^{k + 1}} &\!=\! {\rm{arg}}\mathop {{\rm{min}}}\limits_{{Z}} \left\{ {\lambda {{\left\| {Z} \right\|}_1} \!+\! \left( {\rho /2} \right)\left\| {{{X}^{k \!+\! 1}} \!-\! {Z} \!+\! {{U}^k}} \right\|_{\rm{F}}^2} \right\}\\ & = {S_{\lambda /\rho }}\left( {{{X}^{k + 1}} + {{U}^k}} \right)\\[-12pt] \end{split} $ $

(20)

其中，${S_{\lambda /\rho }}\left( x \right)$表示$x$的软阈值(soft thresholding)算子^[12]，$x$为${X}$中的任意元素，当${X}$属于实数域时，有

$$\begin{split}{{\rm{S}}_{\lambda /\rho }}\left( {{{X}^{k + 1}} + {{U}^k}} \right) =& {\rm{sign}}\left( {{{X}^{k + 1}} + {{U}^k}} \right){\kern 1pt} \\ &\cdot \max \left( {\left| {{{X}^{k + 1}} + {{U}^k}} \right| - \lambda /\rho ,0} \right)\end{split}$ $

(21)

其中，${\rm{sign}}\left( x \right)$与$\max \left( x \right)$分别代表$x$的符号函数与求最大值函数，实数${x_{\rm{R}}}$的软阈值算子${S_{\lambda /\rho }}\left( {{x_{\rm{R}}}} \right)$示意图如图2所示。图2中，实数软阈值可以看作一个分段线性滤波器，对$\left| x \right| \le \left( {\lambda /\rho } \right)$的较小值进行归零操作，当$\left| x \right|$>$\left( {\lambda /\rho } \right)$时，对$x$进行线性缩小。本文针对SAR回波的复数特性，改进现有实数软阈值算子，将复数矩阵${Z}$中的任意元素$z$拆解为实数部分${\rm{RE}}\left( z \right)$和虚数部分${\rm{IM}}\left( z \right)$，并对${X,U}$进行同样复数分解，通过式(20)中的上式对${\rm{RE}}\left( z \right)$和${\rm{IM}}\left( z \right)$分别求导，得到复数软阈值算子，如

图 2  实数软阈值算子示意图

下载: 全尺寸图片 幻灯片

$$\begin{split}{{Z}^{k + 1}} =& {{\rm{S}}_{\lambda /\rho }}\left( {{{X}^{k + 1}} + {{U}^k}} \right)\\ {\rm{ = }}&\left( {{{X}^{k + 1}} + {{U}^k}} \right)/\left| {{{X}^{k + 1}} + {{U}^k}} \right| \\ &\cdot \max (\left| {{{X}^{k + 1}} + {{U}^k}} \right| - \lambda ,0)\end{split}$ $

(22)

由式(22)可知，复数软阈值与实数软阈值一样具有保相性。此处为了表示方便，绘制复数${x_{\rm{C}}}$软阈值算子${S_{\lambda /\rho }}\left( {{x_{\rm{C}}}} \right)$如图3。

图 3  复数软阈值算子示意图

下载: 全尺寸图片 幻灯片

图3中，复数$x$的幅度为$\left| x \right|$，幅角为$\angle x \in \left[ {\left. {0,2{\rm{\pi }}} \right)} \right.$，对应的复数软阈值${{{S}}_{{\lambda /\rho }}}\left( x \right)$的幅角也为$\angle x$，所以由$\left| x \right|$和$\angle x$足以表示所有复数$x$的取值范围, 在平行于纸面的笛卡尔坐标系上画出$\left| x \right|$和${{{S}}_{{\lambda /\rho }}}\left( x \right)$的对应关系也并不困难，然后将坐标系以${{{S}}_{{\lambda /\rho }}}\left( x \right)$轴沿纸面纵旋转整个幅角$\angle x$的取值范围，便描绘出了复数域中$x$与${S_{{\lambda /\rho }}}\left( x \right)$的一一对应关系，实数软阈值也可以看做$\angle x$=$0$时复数软阈值的一个特例。

由于复数软阈值比实数软阈值更具有普适性，尤其适用于复数SAR回波数据的稀疏特征增强处理，基于此，本文提出利用复数软阈值的C-ADMM算法。

步骤 3　对偶变量更新 ${X}$和${Z}$的更新是对其联合取值最小化的过程，在C-ADMM算法中，${U}$的迭代也就对${X}$和${Z}$联合对偶变量的更新过程，即

$${{U}^{k + 1}} = {{U}^k} + {{X}^{k + 1}} - {{Z}^{k + 1}}$ $

(23)

式(19)、式(22)和式(23)给出了C-ADMM更新过程的解析形式，利用高斯-赛德尔思想循环迭代，直到满足停止准则，即达到设定的迭代次数或者精度${\left\| {{{R}^k}} \right\|_{\rm{F}}} \le \varepsilon$，此处${{R}^k} = {{X}^k} - {{Z}^k}$表示主残差，$\varepsilon$为依据经验设定的精度阈值，当满足停止准则时算法终止。

综上所述，可提出基于C-ADMM的SAR稀疏特征增强成像算法流程，如表1所示。

表 1  C-ADMM稀疏特征增强算法流程

(1) 初始化，输入SAR原始数据；

(2) 信号预处理，得到通用信号模型$S\left( {\hat r,t} \right)$或$S\left( {\hat r,{t'}} \right)$；

(3) 设定初值${{X}^0} = {{Z}^0} = {{U}^0} = 0$，构造字典${A}$=${{A}_0}$或${A}\left( {{\gamma _d}} \right)$；

(4) 设定迭代次数与目标精度，若停止准则不满足，进行循环；

(5) 更新目标图像 　 ${{X}^{k + 1}} = {\left( {{{A}^{\rm{H}}}{A} + \rho {I}} \right)^{ - 1}}\left\{ {{{A}^{\rm{H}}}{Y} + \rho \left( {{{Z}^k} - {{U}^k}} \right)} \right\}$;

(6) 更新软阈值${{Z}^{k + 1}} = {S_{\lambda /\rho }}\left( {{{X}^{k + 1}} + {{U}^k}} \right)$;

(7) 更新对偶变量${{U}^{k + 1}} = {{U}^k} + {{X}^{k + 1}} - {{Z}^{k + 1}}$;

(8) 若不满足停止准则，继续步骤5—步骤7，若满足停止准则，跳 出循环；

(9) 输出稀疏特征增强后的图像。

下载: 导出CSV 

| 显示表格

本文在C-ADMM稀疏特征增强算法中，引入增广拉格朗日乘子，进而在每次迭代中应用其作为可行的边界约束，提升了算法的稳健性；通过对偶分解与交替方向迭代，提高了算法的运算效率；利用复数软阈值，增加了针对SAR数据的实用性。

4.   实验验证

为了验证所提的C-ADMM算法与信号模型的有效性和通用性，本节利用复数仿真实验对C-ADMM在不同信噪比和降采样情况下的稀疏恢复性能进行验证，然后使用SAR, SAR-GMTIm和ISAR 3种模式的实测数据进行实验，并与BCS和CVX成像结果进行对比，从而验证C-ADMM稀疏特征增强算法的稳健性和高效性。

4.1   仿真实验

为了定量直观的衡量C-ADMM算法的稀疏恢复性能，本节实验借鉴了Donoho等人^[18]提出的相变图(Phase Transition Diagram, PTD)进行验证。PTD的横坐标为降采样率，纵坐标为稀疏度，图中的每个点代表了长度为1024的原始复数序列在对应降采样率和稀疏度的情况下，通过100组蒙特卡洛实验恢复的成功概率，图中不同颜色表示了平均恢复的成功概率，恢复成功与失败主要通过归一化均方误差(Normalized Mean Square Error, NMSE)进行表征，其表达式为

$${\rm{NMSE}}\; = \;\left\| {\left( {{{X}^k} - {X}} \right)/{X}} \right\|_{\rm{F}}^2$ $

(24)

当NMSE< 0.0001时定义为恢复成功，实验分别对无噪声信号和信噪比为10 dB的噪声信号进行恢复，结果如图4。

图 4  C-ADMM在无噪声和10 db信噪比情况下的相变图

下载: 全尺寸图片 幻灯片

从图4(a)中可以看出，当无噪声时，C-ADMM在较低稀疏度和降采样率的情况下，都能完好恢复原始信号，因此能恢复的区域面积饱满，在稀疏度较大降采样率较低的情况下，即采样率太低以至于无法在采样中保留所有的序列信息，因此难以恢复信号，对应的区域恢复概率接近于0，换言之，原始信号并不具有稀疏性。图4(b)中当加入噪声使信噪比为10 dB时，可以看出能恢复的区域变小，说明加入的噪声影响了部分信号以至于难以恢复，但是当稀疏度较小时仍有大面积的能恢复区域，说明在信号足够稀疏的条件下，C-ADMM仍然具有优异的恢复性能。

4.2   SAR成像模式

传统SAR稀疏特征增强算法运算效率较低^[19,20]，难以适应实际SAR成像应用中较大的成像场景区域和回波数据量。本节应用所提C-ADMM稀疏特征增强算法进行实验验证，实验选用美国Sandia实验室公布的实测SAR复图像数据。该数据应用Ka波段雷达获取，成像分辨率达0.1 m。

图5为应用距离多普勒(Range Doppler, RD)传统算法直接成像与本文中C-ADMM算法成像结果对比，用于参考的RD成像结果背景噪声明显，影响了地面虚线状场景的成像对比度，而C-ADMM成像结果中的强散射点明显被增强，背景噪声受到抑制，地物特征轮廓清晰，虚线部分易于辨别，验证了C-ADMM算法对SAR成像模式稀疏特征增强效果。

图 5  SAR大图C-ADMM稀疏特征增强成像结果

下载: 全尺寸图片 幻灯片

为了进一步评估C-ADMM算法的运算效率和细节表现，图6中实验截取图5部分数据使用C-ADMM与BCS,CVX算法进行成像。如图6(b)所示。应用C-ADMM算法进行稀疏特征增强的细节效果最佳，相比图6(a)用于参考的RD成像结果，背景噪声受到了抑制，稀疏特征得以增强，BCS和CVX算法虽然起到了一定程度的稀疏特征增强效果，但都残留了部分噪声，影响观测结果，C-ADMM成像结果中的虚线目标边界清晰，并且运算效率最高，运算时间$t$仅需5.17 s，相比BCS (t=1481.10 s)和CVX (t=23193.92 s)算法的运算效率提高了2～3个数量级(至少100倍以上)。

图 6  SAR局部C-ADMM稀疏特征增强成像结果

下载: 全尺寸图片 幻灯片

但是由于SAR模式成像区域大，散射特征复杂，不同场景的散射点强度与数量不尽相同，仅调节正则化参数只会抑制弱散射点，保留强散射点，对散射特性不均匀的成像目标很难保证其内部连续性，如图6中C-ADMM在增强稀疏特征的同时，也对线目标的连续性产生了影响，因此在增强稀疏特征的同时，如何有效地保留并增强其他特征，笔者正在探索。

4.3   SAR-GMTIm成像模式

为了验证所提C-ADMM对于地面运动目标的成像性能，本组实验应用美国空军实验室公开的实测Gotcha雷达数据。该数据采用三通道SAR体制，X波段，因此可以初步实现针对地面杂波的抑制效果，成像最高分辨率可达0.3 m。图7(a)所示为两个地面运动目标直接SAR成像结果，由于目标运动，造成成像结果中目标出现严重散焦。本文采用LVD时频表示并构造字典${{A}}$后，分别应用C-ADMM, BCS和CVX进行稀疏特征增强结果如图7(b)—图7(d)所示。对比实验结果可见，C-ADMM和BCS的稀疏特征增强效果较好，目标聚焦性能和图像对比度提升明显。并且C-ADMM计算效率更高(t=5.88 s)，运算耗时明显优于BCS(t=672.12 s)和CVX (t=8657.45 s)。实验结果证明，本文所提C-ADMM稀疏特征增强算法精度完全可以比拟高精度的贝叶斯类算法，在运算效率上相比传统方法至少提升2～3个数量级。

图 7  GMTI 模式C-ADMM稀疏特征增强成像结果

下载: 全尺寸图片 幻灯片

4.4   ISAR成像模式

为了验证信号模型的通用性与C-ADMM算法优势，本组实验使用通过ISAR成像的Yak32的实测数据，ISAR模式一般杂波较少，并且，为了进一步验证C-ADMM的稳健性，本文对数据进行了模拟降采样和加入不同噪声，并用3种方法进行了信号稀疏特征增强，结果如图8与图9。

图 8  ISAR模式C-ADMM不同降采样成像结果

下载: 全尺寸图片 幻灯片

图 9  ISAR模式C-ADMM不同信噪比成像结果

下载: 全尺寸图片 幻灯片

图8在实验数据中均加入了5 dB噪声，从上到下依次为未降采样后的图像、降采样率为0.5和0.25的图像，从左到右依次为RD, C-ADMM, BCS, CVX处理的结果。从RD成像结果图8(a)、图8(e)、图8(i)看出，降采样率越低，成像时产生的噪声越大，恢复的图像效果越差，飞机轮廓难以识别。其它3种算法成像结果都对噪声进行了抑制，其中在同样降采样情况下(降采样0.5)，C-ADMM和BCS的稀疏特征增强效果最好，背景噪声最为干净，目标飞机轮廓清晰，其中，在处理时间t上C-ADMM更为高效(t=0.41 s)，比BCS (t=93.88 s)和CVX (t=528.55 s)快了至少2个数量级。

图9对实验数据均进行了0.5倍降采样，从上到下的信噪比情况依次为5 dB, 0 dB, –2 dB，在参考的 RD算法结果中，随着加入噪声越来越大，成像目标难以辨识，直至淹没在噪声中。从图中可以看出，C-ADMM等3种都产生了稀疏特征增强效果，当噪声较大时，BCS在部分距离单元上的噪声没有抑制效果，而CVX成像结果丢失了目标的部分轮廓，C-ADMM的稀疏效果最好，完好的保留了目标的轮廓，并抑制了背景噪声，在处理时间t上，C-ADMM相比BCS和CVX快了90倍以上。

在以上3组实验中，由于SAR与GMTI模式强散射点，成像场景噪声情况复杂，可以使用C-ADMM算法在成像处理的同时，对成像目标的稀疏特征进行增强，而ISAR模式中成像目标相对于背景本身就具有强稀疏性，因此可以使用C-ADMM算法进行稀疏信号恢复，多种噪声与降采样情况的恢复结果都直观地反映了C-ADMM算法的稳健性。在成像过程中，CVX算法利用递归最小二乘思想迭代求解，而BCS每次迭代必须进行求逆运算，这导致了CVX和BCS方法的求解效率并不高，C-ADMM结合对偶思想并引入拉格朗日乘子，求解效率更高，相比BCS和CVX，在处理效率上节省了近百倍的时间，说明C-ADMM在SAR稀疏特征增强中更为高效。

5.   结束语

本文针对SAR, GMTIm和ISAR 3种成像模式，推导获得了通用的回波信号模型，同时提出了面向复数回波数据的C-ADMM稀疏特征增强算法，可实现目标稀疏特征的稳健高效恢复。本文实验应用多组数据进行验证，同时对比传统CVX和BCS稀疏特征增强方法，证明了本文所提C-ADMM算法的有效性和优越性。然而，在实验过程中，本文发现单一增强稀疏特征对目标成像精度的提升帮助有限，有必要研究适用于更多特征的增强方法，这将是后续工作的重点。