1 引言

逆合成孔径雷达(Inverse Synthetic Aperture Radar,ISAR)可在远距离、全天候条件下对运动目标进行成像，因而在战略预警、目标识别以及空中交通管制等军事和民用领域发挥着重要的作用，并得到了广泛的应用^[1]。ISAR能够同时获得目标距离向和方位向的高分辨，从而可使雷达目标呈现2维图像。其中，距离向高分辨率取决于信号的带宽，带宽越大，距离分辨率越高；方位向高分辨率则通过提高方位向的积累时间(目标沿雷达视线方向的转角)获得。然而，各种不利条件对传统ISAR成像带来了挑战，比如回波被强噪声污染或由于敌方干扰需要人为地删除部分数据造成回波缺损时，如何保证高分辨率及高质量成像效果。

近年来压缩感知(Compressive Sensing,CS)理论^[2]的出现为上述问题的解决提供了一种新思路。由于ISAR成像满足稀疏性机理^[3]，可以将CS理论应用于ISAR成像。自Baraniuk将CS理论应用于雷达成像以来^[3]，国内外众多学者开展了基于CS理论3个关键环节(雷达回波稀疏表示、量测矩阵构造和重构算法设计)的ISAR成像研究。文献^[4]针对ISAR 信号在传统傅里叶基下稀疏表示时出现的能量泄露和频谱展宽的问题，提出了一种改进的傅里叶基，并通过仿真表明改进稀疏基的成像效果优于传统傅里叶基。文献^[5]研究了基于ISAR信号的参数化稀疏表示模型，利用求解该模型的稀疏性算法既能补偿目标旋转速度又能进行ISAR成像，从旋转速度估计的准确度、ISAR成像的效果和计算的复杂度3个方面比较了4种参数化稀疏重构算法，给出了每种算法的优点和缺点。文献^[6]考虑到散射点之间的相关性，联合局部稀疏性约束和非局部全变分提出一种ISAR成像方法，该方法能够有效地压制噪声及去除虚假点，缺点是算法运算量较大。文献^[7] 基于CS提出了一种2维联合超分辨ISAR成像算法， 由于该方法有效利用了2维耦合信息，因此提高了2维分辨率。基于CS理论的ISAR成像方法已涌现出很多成果，但是，目前CS理论在ISAR成像中的应用大多数都是基于单量测向量(Single Measurement Vector,SMV)模型，即将ISAR回波数据按照距离单元划分，对每一距离单元重构后组合得到最终的成像结果。由于SMV模型下，基于CS理论的ISAR 成像重构算法效率较低及实时性较差，使得CS理论距离实际应用还有较大差距；而近年来出现的基于多量测向量(Multiple Measurement Vectors,MMV) 模型^[8]的CS理论可在一定程度上减少重构时间，提高重构质量。但是MMV模型需要信号满足联合稀疏的假设，文献^[9]将相位自聚焦模型等效为MMV 模型，基于多任务贝叶斯模型提出了一种相位自聚焦ISAR成像算法，取得了较好的成像效果。文献^[10] 将SAR回波在距离向上等效为MMV模型，实现了快速SAR成像。而这在ISAR信号模型中一般是不满足的。因此研究符合ISAR回波信号特点的MMV 模型，提高ISAR成像效率是一个值得研究的问题。

本文针对这一问题，首先依据多量测向量的ISAR成像模型，研究了基于线性Bregman迭代(Linearized Bregman Iteration,LBI)的一类多量测向量成像算法。全文安排如下：第2节介绍基于多量测向量的ISAR成像模型；并针对该模型研究了一类基于LBI的多量测向量成像算法，给出了此类算法的整体构架和应用条件，然后给出了收敛性和抗噪性分析；第3节通过仿真结果分析了算法的良好性能，分析了正则化参数对算法性能的影响，并通过实测数据的实验结果验证了所提算法的有效性及优势；最后给出总结。

2 基于线性Bregman迭代类的多量测向量成像算法 2.1 基于多量测向量的ISAR成像模型

本文的研究重点是ISAR方位向成像，假设已对回波数据完成了脉冲压缩、平动补偿^[11]。假设目标在X-Y 2维平面内绕自身重心原点O以角速度w逆时针匀速转动，如图 1所示。

设雷达距目标坐标系原点O的距离为R_o，目标上某一散射点P在初始时刻位于目标直角坐标系的坐标为(x_P,y_P)，散射点P到目标坐标系原点O的距离为${R_{\rm{P}}} = \sqrt {x_{\rm{P}}^2 + y_{\rm{P}}^2} $,P点初始相位为 θ，观测时间为T_a。t_P时刻P点到雷达的瞬时距离为：

$\begin{align} &R({{t}_{\text{P}}})=\sqrt{R_{\text{o}}^{2}+R_{\text{P}}^{2}+2\cdot {{R}_{\text{o}}}\cdot {{R}_{\text{P}}}\cdot \sin (\theta +\omega {{t}_{\text{P}}})} \\ &\approx {{R}_{\text{o}}}+{{x}_{\text{P}}}\omega {{t}_{\text{P}}}+{{y}_{\text{P}}} \\ \end{align}$

(1)

假设雷达发射线性调频信号

$s(t)=\text{rect}\left( \frac{{{t}_{\text{f}}}}{{{T}_{\text{p}}}} \right)\cdot \text{exp}\left( \text{j}2\pi \left( {{f}_{\text{c}}}t+\frac{1}{2}{{k}_{\text{r}}}t_{\text{f}}^{2} \right) \right)$

(2)

其中， μ>0为矩形函数，f_c为载频，T_p为脉宽，k_r为调频率，t_f=t–nT为距离向快时间，t_n=nT为方位向慢时间，T为脉冲重复周期。则P点的回波经距离向脉压后为：

$\begin{align} &{{s}_{\text{r}}}({{t}_{\text{f}}},{{t}_{n}})=\sigma \text{rect}\left( \frac{{{t}_{n}}}{{{T}_{\text{a}}}} \right)\text{sinc}\left[ {{T}_{\text{p}}}{{k}_{\text{r}}}\left( {{t}_{\text{f}}}-\frac{2\left( {{R}_{\text{o}}}+{{y}_{\text{P}}} \right)}{c} \right) \right] \\ &\cdot \exp \left[ -\text{j}4\pi \frac{\left( {{R}_{\text{o}}}+{{y}_{\text{P}}} \right)}{\lambda } \right]\exp \left( -\text{j}2\pi {{f}_{d}}{{t}_{n}} \right) \\ \end{align}$

(3)

其中， σ为P点的散射系数，c为光速， λ为波长，多普勒频率f_d=2 wx_P/ λ。此时即实现了距离向分辨。将经过距离对准后的回波按距离单元排列，假设共取D个距离单元，若第d个距离单元内有K个强散射点，则此距离单元的回波信号可表示为：

${{s}_{d}}\left( {{t}_{n}} \right)=\sum\limits_{k=1}^{K}{{{\sigma }_{k}}}\cdot \text{rect}\left[ \frac{{{t}_{n}}}{{{T}_{\text{a}}}} \right]\exp \left( -\text{j}2\pi {{f}_{k}}{{t}_{n}} \right)$

(4)

其中，f_k为第k个散射点的多普勒频率。由式(4)可知，与整个多普勒域相比，目标所占的多普勒成分仅有K个，满足稀疏性，因此可以构建基于CS的ISAR模型。

定义慢时间间隔为 Δt以及多普勒频率分辨率为 Δf_d，假设脉冲重复频率为f_r，则慢时间序列t_n可以表示为${{t}_{n}}={{\left[ 1:{{N}_{\text{a}}} \right]}^{\text{T}}}\Delta t,\Delta t=1/{{f}_{\text{r}}}$，其中N_a=T/ Δt表示脉冲数。假设多普勒域范围[–f_r/2,f_r/2]，则多普勒序列为${{f}_{d}}=\left[ 1:N \right]\Delta {{f}_{d}}-{{f}_{\text{r}}}/2$,N=f_r/ Δf_d，其中N为多普勒单元数。位于各多普勒单元散射点的强度用向量x_d表示x_d=[x_d(1) x_d(2)  ··· x_d(N)]^T，方位向稀疏基可构造为 Ψ=exp(–j2 πt_nf_d)， 故可将第d个距离单元的方位成像问题等价为如下重构问题：

$\begin{align} &\min {{\left\| {{x}_{d}} \right\|}_{0}}, \text{s}.\text{t}. {{y}_{d}}\text{=}\Phi {{s}_{d}}\left( {{t}_{n}} \right)=\Phi \Psi {{x}_{d}}=\Theta {{x}_{d}}, \\ &d=1,2,\cdots ,D \\ \end{align}$

(5)

其中， Φ为M×N_a的量测矩阵，M为量测数， Θ= Φ Ψ为M×N的感知矩阵。文献^[2]证明了当感知矩阵 Θ满足RIP条件时，信号可以从量测值y_d中精确重构出来，此时可以利用CS重构算法求解，对D个距离单元逐一重构，然后组合即可获得ISAR 图像。由于对每个距离单元都需要重构，将会产生大量的冗余计算，增加不必要的计算负担^[12]。

为减少冗余计算，缩短成像时间，本文研究基于多量测向量的ISAR成像模型，即

$Y=\Theta X$

(6)

其中，Y=[y₁ y₂ ··· y_D]为回波数据矩阵，X=[x₁ x₂ ··· x_D]为待重构的散射点矩阵。若求得解X，即可得到ISAR图像。与传统MMV模型不同，这里不同距离单元散射点的稀疏度和支撑集都不需要相同。令k_d表示x_d的稀疏度，向量K=[k₁ k₂ ··· k_D]表示X的稀疏度，称为X的稀疏度向量，即X为K稀疏的，此时X不再是联合稀疏的。首先利用X的稀疏特性约束，将求解式(6)的稀疏解问题描述为如下优化问题：

$\underset{X\in {{C}^{N\times D}}}{\mathop{\min }}\,{{\left\| X \right\|}_{0,q}}, \text{s}.\text{t}. Y=\Theta X$

(7)

其中，${{\left\| X \right\|}_{0,q}}=\left| \text{supp}\left( X \right) \right|$,supp(X)为X的支撑集，即非零元素的个数。

式(7)是NP难的问题，在感知矩阵 Θ满足RIP条件时，可以将式(7)转化为以下凸优化问题：

$\underset{X\in {{C}^{N\times D}}}{\mathop{\min }}\,{{\left\| X \right\|}_{1}}, \text{s}.\text{t}. Y=\Theta X$

(8)

其中，矩阵的1-范数定义为${{\left\| X \right\|}_{1}}=\sum\limits_{i=1}^{N}{\sum\limits_{j=1}^{N}{\left| {{X}_{i,j}} \right|}}$。

2.2 算法原理

有关LBI理论的详细内容可参考文献^[13]，文献^[13]首先提出Bregman迭代方法，并提出LBI。此后，文献^[14-16]深入研究了LBI的性能，取得了一系列的研究成果。

本节首先根据多量测向量ISAR成像模型，然后利用LBI理论的基本原理，推导出了多量测向量线性Bregman算法(Multiple measurement vectors Linearized Bregman,MLB)算法。针对MLB算法存 在 “停滞现象 ”，收敛速度慢的问题，研究了多量测向 量快速线性Bregman算法(Multiple measurement vectors Fast Linearized Bregman,MFLB)算法、多量测向量改进线性Bregman算法(Multiple measurement vectors Improved Linearized Bregman,MILB)算法和多量测向量快速改进型线性Bregman算法(Multiple measurement vectors Fast Improved Linearized Bregman,MFILB)算法以加快收敛速度，减少迭代次数，缩短成像时间。为便于理解，本文首先给出这4种算法的相互关系及适用条件，如表 1所示，然后进一步分析各类算法特征。从表中可以看出，基于LBI理论得到MMV-ISAR模型的成像算法，即得到MLB算法后，首先判断感知矩阵的条件数是否为1，若为1，则继续判断是否存在停滞，若不存在，则此时使用MLB算法进行成像，若存在停滞，则此时使用MFLB算法进行成像；若感知矩阵条件数不为1，则继续判断是否存在停滞，若不存在，则此时使用MILB算法进行成像，若存在停滞，则此时使用MFILB算法进行成像。

为从含有噪声回波数据矩阵中恢复X，放松式(8)的约束项，利用正则化参数 m控制稀疏度和误差，转化为如下正则化形式：

$\widehat{X}\text{=}\underset{X\in {{C}^{N\times D}}}{\mathop{\text{argmin}}}\,\mu {{\left\| X \right\|}_{1}}\text{+}\frac{1}{2}\left\| \Theta X-Y \right\|_{\text{F}}^{2}$

(9)

其中，${\left\| { \cdot } \right\|_{\rm{F}}}$表示矩阵的Frobenius范数，简称F范数，定义如下：

${{\left\| A \right\|}_{\text{F}}}=\sqrt{\left( \sum\limits_{i,j=1}^{N}{{{\left| {{a}_{ij}} \right|}^{2}}} \right)}=\sqrt{\text{tr}\left( {{A}^{\text{H}}}A \right)}$

(10)

其中，tr(·)为矩阵的迹，即对角线元素之和。

2.2.1 MLB 算法

本文成像模型与文献^[17]中的求解模型是一致的，我们在文献^[17]中详细推导了算法原理，下面给出MLB算法的迭代格式：

$\left\{ \begin{array}{*{35}{l}} {{V}^{k+1}}={{V}^{k}}+{{\Theta }^{\text{H}}}(Y-\Theta {{X}^{k}}) \\ {{X}^{k+1}}\text{=}\delta \text{sof}{{\text{t}}_{\mu }}({{V}^{k+1}}) \\ \end{array} \right.$

(11)

其中，soft(·)为复数矩阵条件的软阈值算子，定义如下：

$\begin{align} &\text{sof}{{\text{t}}_{\mu }}\left( V \right)=\frac{\max \left\{ \left| {{V}_{ij}} \right|-\mu ,0 \right\}}{\max \left\{ \left| {{V}_{ij}} \right|-\mu ,0 \right\}+\mu }{{V}_{ij}} \\ &=\left\{ \begin{matrix} \frac{{{V}_{ij}}}{\left| {{V}_{ij}} \right|}\left( \left| {{V}_{ij}} \right|-\mu \right),&\left| {{V}_{ij}} \right|>\mu \\ 0,&\left| {{V}_{ij}} \right|\le \mu \\ \end{matrix} \right. \\ \end{align}$

(12)

其中，V_ij为复数矩阵V中第i行第j列的元素，|V_ij|表示V_ij的模。从软阈值算子soft的定义中可以看出，MLB算法的软阈值算子从点收缩拓展为面收缩，平面阈值比原始的线阈值具有更好的压制噪声效果，因为平面阈值考虑到了全局数据，而线阈值仅考虑到当前一列的数据，基于实测数据的成像实验将验证这一结论。

此外，由于MLB算法只是最基本的迭代格式，一方面，存在 “停滞 ”现象，即几次迭代时，算法输出保持不变；另一方面，MLB算法在感知条件数较大时，收敛速度将会较慢。为解决这两方面的问题，研究3种算法，下面具体分析。

2.2.2 MFLB 算法

文献^[14]通过估计迭代步长的方 法减少了LB算法的 “停滞 ”现象带来的迭代次数变大的影响，下面利用这一思想对MLB算法实现快速迭代。在停滞时间内，假设发生m步停滞，第k次迭代时停滞开始，第k+m次迭代时停滞终止，第k+m+1次迭代时V^k+m+1突破了闭区间[– μ, μ]，产生一个新的X^k+m+1(X^k+m+1≠X^*)，因为X^k+m最终能够收敛到一个极限值X^*[14]，因此有X^k${{X}^{k}}\approx {{X}^{k+1}}\approx {{X}^{k+2}}\approx \cdots {{X}^{k+m}}\approx {{X}^{*}}$，所以V的增量 ΔV=Θ^H(Y– ΘX^k)可认为是固定的，那么停滞期间的迭代过程可表示为：

$\left\{ \begin{array}{*{35}{l}} {{X}^{k+j}}\equiv {{X}^{k}} \\ {{V}^{k+j}}\text{=}{{V}^{k}}+j\cdot \Delta V,\ j=1,\cdots ,m \\ \end{array} \right.\ \ $

(13)

在停滞期间内，每次迭代V^k仅以固定的增量 ΔV进行积累，第j次迭代的积累量为j· ΔV,V^k+j不足以突破闭区间[–μ, μ]，第k+m+1次迭代的V^k+m+1才能突破闭区间[–μ, μ]。因此，要想终止停滞，需要将V^k的积累量增大为(k+m+1)· ΔV，才能使V^k+m+1能够突破闭区间[–μ, μ]的限制。

为计算停滞步长，首先对X^k,V^k,X^k+j,V^k+j 和 ΔV进行矩阵向量化处理，即${{\widetilde{X}}^{k}}\text{=vec}\left( {{X}^{k}} \right),{{\widetilde{X}}^{k+j}}=\text{vec}\left( {{X}^{k+j}} \right),{{\widetilde{V}}^{k+j}}=\text{vec}\left( {{V}^{k+j}} \right)$和$\Delta \tilde{V}=\text{vec}\left( \Delta V \right)$。定义I₀为${{\widetilde{X}}^{k}}$零元素的索引集，${I_1} = {\bar I_0}$为${{\widetilde{X}}^{k}}$非零元素的支撑集，此时式(13) 可改写为如下分段形式：

$\left\{ \begin{matrix} \widetilde{X}_{i}^{k+j}\equiv \widetilde{X}_{i}^{k},&\forall i \\ \widetilde{V}_{i}^{k+j}\text{=}\widetilde{V}_{i}^{k}+j\cdot \Delta {{\widetilde{V}}_{i}},\ &i\in {{I}_{0}} \\ \widetilde{V}_{i}^{k+j}\equiv \widetilde{V}_{i}^{k},\ &i\in {{I}_{1}} \\ \end{matrix} \right.$

(14)

当且仅当I₀中${{\widetilde{V}}^{k}}$的元素突破闭区间[–μ, μ]的限制时，${{\widetilde{X}}^{k}}$才会产生一个新的非零元素，从而消除停滞现象。当$i\in {{I}_{0}},\widetilde{V}_{i}^{k}\in \left[ -\mu ,\mu \right]$,时，可以利用式(15)估计$\widetilde{V}_{i}^{k}$突破限制需要的积累步数，即

$s=\text{min}\left\{ {{s}_{i}}\left| {{s}_{i}}=\left\lceil \frac{\mu \cdot \text{sign}\left( \Delta \widetilde{V}{{}_{i}} \right)-\widetilde{V}_{i}^{k}}{\Delta {{\widetilde{V}}_{i}}} \right\rceil ,\forall i\in {{I}_{0}} \right| \right.$

(15)

其中，$\left\lceil \centerdot \right\rceil $为取整符号，s为需要的积累步数，也就是 停滞的长度。得到s后就可以利用式(16)终止停滞。

$\left\{ \begin{array}{*{35}{l}} {{X}^{k+s}}\equiv {{X}^{k}} \\ {{V}^{k+s}}\text{=}{{V}^{k}}+s\cdot \Delta V \\ \end{array} \right.$

(16)

因此，当X^k在两步迭代保持不变时，认为其处于迭代停滞状态，此时可通过增加V^k的变化量以突破[–μ, μ]，从而使X^k加速到停滞的临界点，以此减少积累时间，加快算法运行速度。

注意到，MFLB算法虽然通过估计迭代步长加快了算法收敛，但依然没有解决感知条件数较大时，收敛速度将较慢的问题。

2.2.3 MILB 算法

文献^[15]指出，感知矩阵条件数越小，则LB算法收敛越快。因此，需要对感知矩阵的条件数进行优化以减少迭代次数，具体的优化方法在文献^[17]中进行了详细的介绍。下面介绍关键的步骤。

首先在式(6)两边乘以( Θ Θ^H)^–1/2，得到

${{(\Theta {{\Theta }^{\text{H}}})}^{-1/2}}Y={{(\Theta {{\Theta }^{\text{H}}})}^{-1/2}}\Theta X$

(17)

令${{Y}_{1}}={{(\Theta {{\Theta }^{\text{H}}})}^{-1/2}}Y,{{\Theta }_{1}}={{(\Theta {{\Theta }^{\text{H}}})}^{-1/2}}\Theta $，则有

${{Y}_{1}}={{\Theta }_{1}}X$

(18)

对式(18)采用MLB算法求解，将Y₁, Q₁代入，可得MIBL算法迭代公式

$\left\{ \begin{array}{*{35}{l}} {{V}^{k+1}}={{V}^{k}}+{{\Theta }^{\text{H}}}{{(\Theta {{\Theta }^{\text{H}}})}^{-1}}(Y-\Theta {{X}^{k}}) \\ {{X}^{k+1}}=\delta \text{sof}{{\text{t}}_{\mu }}({{V}^{k+1}}) \\ \end{array} \right.$

(19)

其中，${{\Theta }^{\text{H}}}{{(\Theta {{\Theta }^{\text{H}}})}^{-1}}$是Θ的右逆，当 Θ不满秩时，用 Θ的Moore-Penrose逆 Θ⁺代替右逆。

由式(18)可知此时感知矩阵为 Θ₁，计算 Θ₁条件数

${{\Theta }_{1}}\Theta _{1}^{\text{H}}=({{(\Theta {{\Theta }^{\text{H}}})}^{-1/2}}\Theta ){{({{(\Theta {{\Theta }^{\text{H}}})}^{-1/2}}\Theta )}^{\text{H}}}=I$

(20)

$\text{cond}({{\Theta }_{1}})={{[{{\lambda }_{\max }}({{\Theta }_{1}}\Theta _{1}^{\text{H}})/{{\lambda }_{\min }}({{\Theta }_{1}}\Theta _{1}^{\text{H}})]}^{1\text{/2}}}=1$

(21)

可以看出MILB算法对感知矩阵的条件数进行了优化，达到了加快收敛速度的目的。

2.2.4 MFILB 算法

MFLB算法和MILB算法在一定程度上能够减少迭代次数，但这两种算法都存在缺陷：一是MFLB算法虽然通过估计迭代步长减少了迭代次数，但是该算法的感知矩阵条件数不是最优的，所以迭代次数不是最快的；二是MILB算法虽然将感知矩阵条件数优化到最小，但该算法依然存在 “停滞 ”现象。因此本节从这两方面考虑，首先对感知矩阵的条件数进行优化，从整体上加快了收敛的速度，然后对停滞步长进行估计，以进一步提高收敛速度，具体内容可参考文献^[18]。

优化感知矩阵条件数的方法与MILB算法相同，估计迭代步长的方法与MFLB的思想大体一致， 需要注意的是MFILB算法中 ΔV= Θ⁺(Y– QX^k)。 具体过程与MFLB算法类似，得到s，也就是停滞的长度，然后利用式(23)终止停滞。

$s\text{=min}\left\{ {{s}_{i}}\left| {{s}_{i}}\text{=}\left\lceil \frac{\mu \cdot \text{sign}\left( \Delta \widetilde{V}{{}_{i}} \right)-\widetilde{V}_{i}^{k}}{\Delta \widetilde{V}{{}_{i}}} \right\rceil ,\ \forall i\in {{I}_{0}} \right. \right\}$

(22)

$\left\{ \begin{array}{*{35}{l}} {{X}^{k+s}}\equiv {{X}^{k}} \\ {{V}^{k+s}}\text{=}{{V}^{k}}+s\cdot \Delta V \\ \end{array} \right.$

(23)

2.3 性能分析 2.3.1 收敛性分析

MLB算法、MFLB算法、MILB算法和MFILB算法是LBI在MMV模型下的拓展，因此分析算法收敛性时，重点分析其是否依然满足收敛性结论。MMV模型下，可将迭代格式改写为如下形式：

$\left\{ \begin{array}{*{35}{l}} \left[ \begin{matrix} V_{1}^{k+1}&\cdots &V_{D}^{k+1} \\ \end{matrix} \right] \\ =\left[ \begin{matrix} V_{1}^{k}&\cdots &V_{D}^{k} \\ \end{matrix} \right]+C(\left[ \begin{matrix} Y{{}_{1}}&\cdots &Y{{}_{D}} \\ \end{matrix} \right] \\ \quad -\Theta \left[ \begin{matrix} X_{1}^{k}&\cdots &X_{D}^{k} \\ \end{matrix} \right]) \\ \left[ \begin{matrix} X_{1}^{k+1}&\cdots &X_{D}^{k+1} \\ \end{matrix} \right] \\ =\delta \text{sof}{{\text{t}}_{\mu }}\left( \left[ \begin{matrix} V_{1}^{k+1}&\cdots &V_{D}^{k+1} \\ \end{matrix} \right] \right) \\ \end{array} \right.$

(24)

其中，当C= Q^H时，迭代式(24)为MLB和MFLB；当C= Q⁺时，迭代式(24)为MILB和MFILB。则式(24)的第j列为：

$\left\{ \begin{array}{*{35}{l}} \left[ V_{j}^{k+1} \right]=\left[ V_{j}^{k} \right]+C\left( \left[ {{Y}_{j}} \right]-\Theta \left[ X_{j}^{k} \right] \right) \\ \left[ X_{j}^{k+1} \right]=\delta \text{sof}{{\text{t}}_{\mu }}(\left[ V_{j}^{k+1} \right]) \\ \end{array} \right.$

(25)

式(25)是SMV模型的LBI算法，每一列满足收敛性结论，因此拓展到MMV后，依然满足收敛性结论。

2.3.2 抗噪性分析

下面对所提4种算法的抗噪性能进行分析，由于3种改进算法只是在收敛速度得到了加快，因此4种算法具有相同的抗噪性能，下面进行统一分析。

下面对抗噪性能进行具体分析。为方便解释抗噪性能，将迭代式(11)改写为如下等价形式：

$\left\{ \begin{array}{*{35}{l}} {{Y}^{k+1}}=Y+{{Y}^{k}}-\Theta {{X}^{k}} \\ {{X}^{k+1}}=\delta \text{sof}{{\text{t}}_{\mu }}({{\Theta }^{\text{H}}}{{Y}^{k+1}}) \\ \end{array} \right.$

(26)

当k=0,Y⁰=0,X⁰=0时，Y¹=Y，将迭代式(26) 中的含噪输入Y¹分解为两部分：Y¹= QX¹+ QB¹，其中X¹可看作为原始纯净信号X的一部分，收缩算子soft可将 Θ^HY¹中的小信号成分过滤掉，因此X¹是过平滑的且不含任何噪声；B¹包含两部分：一是原始纯净信号X中未恢复的信号X-X¹；二是噪声分量Ω，可表示为 ΘB¹= Θ(X-X¹)+ Ω，又有Y¹= ΘX¹+ ΘB¹，所以Y¹= Θ(X-X¹)+ Ω+ ΘX¹。若期望从B¹中恢复出未恢复信号X-X¹，则需要在第2次迭代时将 ΘB¹反馈到原始含噪量测输入Y中，所以第2次迭代新的输入Y²为：

$\begin{align} &{{Y}^{2}}=Y+\Theta {{B}^{1}}=2\Theta {{X}^{1}}+2\Theta {{B}^{1}}-\Theta {{X}^{1}} \\ &=2\Theta \left( \overline{X}-{{X}^{1}} \right)+2\Omega +\Theta {{X}^{1}} \\ \end{align}$

(27)

和第1次迭代含噪量测输入Y¹相比，未恢复信号X-X¹变为两倍，同时Y²包含的噪声分量也变为两倍。由于第2次迭代新的含噪量测输入Y²可分解为Y²= ΘX²+ ΘB²，利用Y²求解X²时，Y²中的信号成分不仅使X²继承了X¹，而且重构了未恢复信号X-X¹的部分信息，因此X²比X¹更逼近X。那么经过k次迭代后，有${{Y}^{2}}=k\Theta \overline{X}+k\Omega -\sum\nolimits_{j}^{k-1}{\Theta }{{X}^{j}}$,X^k能够逼近于X-X¹ ，虽然X^k最终会含有部分噪声，然而通过选择最佳的 μ，可使得${{\left\| Y-\Theta {{X}^{k}} \right\|}_{\text{F}}}$逼近于${{\left\| Y-\Theta \overline{X} \right\|}_{\text{F}}}$，表明算法具备良好的抗噪性能。若已知总的噪声方差∑²，则可以利用式(28)作为噪声条件下的停止准则：

$\left\| Y-\Theta {{X}^{k}} \right\|>\left\| Y-\Theta \overline{X} \right\|={{\sum }^{2}}$

(28)

3 实验与分析 3.1 性能与分析

本文所有实验数据处理都是在操作系统为Windows7的个人计算机上实现的，仿真平台为Matlab R2008b，计算机主要参数如下：处理器为Intel酷睿E7500，主频为2.93 GHz，内存为2 GB。

仿真 1 迭代次数及运算时间比较

本仿真的主要目的是比较LB,FLB,ILB和FILB 4种算法迭代次数及运算时间的关系。参数设置：稀疏信号长度为400，稀疏度为40，量测矩阵为200行400列的随机高斯矩阵，停止准则为${{\left\| \Theta {{X}^{k}}-Y \right\|}_{\text{F}}}/{{\left\| Y \right\|}_{\text{F}}}\le {{10}^{-5}}$，相对重构误差定义为：$\text{error}={{\left\| {{X}^{k}}-X \right\|}_{\text{F}}}/{{\left\| X \right\|}_{\text{F}}}$。仿真结果如图 2 及表 2所示。

由图 2及表 2可以看出，4种算法的相对重构误差随迭代次数的增加，最终都能减小到设定的停止门限。同时MLB算法存在明显的 “停滞现象 ”；MFLB算法通过估计迭代步长的方式有效地缓解了 “停滞现象 ”；MILB算法通过优化感知矩阵条件数的方式有效地减少了迭代次数，加快了收敛；MFILB 算法在优化感知矩阵条件数的基础上估计迭代步长，进一步减少了迭代次数，体现出MFILB算法的优势。

仿真 2 正则化参数对算法性能的影响

本仿真主要验证正则化参数对4种线性Bregman 算法性能的影响。根据文献^[15]的收敛性定理和本文的收敛性分析可知，0< d<1$0<\delta <1/\left\| \Theta {{\Theta }^{H}} \right\|$时MLB和MFLB算法能够保证收敛；0< d<1时MILB和MFILB算法能够保证收敛。因此在仿真中固定MLB和MFLB算法的 d=0.9$\delta \text{ = 0}\text{.9}/\left\| \Theta {{\Theta }^{H}} \right\|$，固定MILB和MFILB算法的 δ=0.9。信号长度N=400,$\mu =k{{\left\| {{\Theta }^{*}}Y \right\|}_{\infty }}$ (MLB和MFLB算法中 Θ^*= Θ^H,MILB和MFILB算法中 Θ^*= Θ⁺),k=[1:3:19]。为综合比较算法性能，本文引入相变图^{[19, 20]}进行衡量，定义采样率 d为量测值和信号长度的比值，即 δ=M/N取值范围是[0.05,0.95]，共取20个等间隔点。将稀疏度和量测值的比值定义为 ρ，即 ρ=K/M，取值范围是[0.05,0.95]，同样取20个等间隔点。每对 δ和 ρ对应一个相对重构误差errer，完成30次独立重复实验，若errer<0.01，则认为此次重构成功，记录成功重构的概率，即( δ, ρ,P)，据此画出相变图，为方便表示，将P=0.5的曲线画出进行比较。图 3～图 6中图(a)分别为MLB,MFLB,MILB和MFILB算法不同 m的相变图，图 3～ 图 6中图(b)分别为对应的平均迭代次数，平均迭代次数的计算方式为：先对重复试验次数取平均，然后再对 ρ取平均。

根据图 3～图 6， m较小时，4种线性Bregman 算法重构性能都较差，随着 m越大，4种算法重构性能都越好，但所需迭代次数也越多，因此，在选择 m时既不能过小( m过小导致重构性能较差)，又不能过大( m过大导致所需迭代次数过多，可行性较差)， 在实际中可依据具体的需求选择不同的正则化参数。同时。可以看出，所需平均迭代次数的关系为：MFILB<MILB<MFLB<MLB，此结论与仿真1也是一致的。

仿真 3 与其他算法的比较

本仿真将4种线性Bregman算法与正交匹配追踪(Orthogonal Matching Pursuit,OMP)算法^[21]、IST^[22]和IHT^[23]算法进行比较。信号模型参数设置与仿真2相同，OMP算法中稀疏度冗余设置为2K； 本文方法设置 m=19$\mu =19{{\left\| {{\Theta }^{*}}Y \right\|}_{\infty }}$，结果如图 7所示。

从图 7可以看出，IHT算法性能优于IST算法；OMP算法性能优于IHT算法；MLB和MFLB算法的性能优于OMP算法；MILB和MFILB算法的性能优于MLB和MFLB算法。

3.2 实测数据成像实验

为进一步验证算法性能，下面利用实测数据进行不同信噪比的ISAR成像实验。部分雷达参数如下：雷达发射信号为LFM信号，带宽为100 MHz， 工作频段为S波段。图 8(a)为512脉冲RD成像结果，也是后续超分辨重建算法的参考模板，超分辨重建算法利用128脉冲，然后进行4倍超分辨成像；图 8(b)为目标支撑区，是通过合适的门限将图 8(a) 划分为目标支撑区和背景支撑区，对比算法得到的超分辨图像若落在目标支撑区则判定为目标，反之则判定为噪声。实验中通过以下3个评价指标对成像质量进行定量分析。

(1) 目标背景比(Target-to-Background Ratio,TBR)。定义为：

$\text{TBR}=10\log 10\left( \frac{\sum\limits_{\left( i,j \right)\in T}{{{\left| {{\widehat{X}}_{i,j}} \right|}^{2}}}}{\sum\limits_{\left( i,j \right)\in B}{{{\left| {{\widehat{X}}_{i,j}} \right|}^{2}}}} \right)$

(29)

其中，T和B分别表示目标支撑区和背景支撑区，$\widehat{X}$表示通过重建算法得到的超分辨图像。TBR表征了重建算法的估计精度和对噪声的抑制能力^[24]。 TBR越大则估计精度越高、对噪声的抑制能力越强。

(2) 噪声能量(Noise energy)。定义为：

$\text{Noise}\ \text{energy}=\sum\limits_{\left( i,j \right)\in B}{{{\left| {{\widehat{X}}_{i,j}} \right|}^{2}}}$

(30)

噪声能量表示落入背景支撑区的虚假散射点和残余噪声能量水平，噪声能量越大，目标识别率越低。

(3) 相关系数(Coherence)。定义为：

$\text{Coherence}=\frac{\left\langle \left| X \right|\cdot \left| \widehat{X} \right| \right\rangle }{\sqrt{{{\left| X \right|}_{2}}}\sqrt{{{\left| \widehat{X} \right|}_{2}}}}$

(31)

其中，$\widehat{X}$为参考模板，也就是图 8(a)所示的512脉冲的RD图像。相关系数表示了超分辨图像与参考模板之间的相似程度，相关系数越大，在基于模板匹配的雷达目标识别中能获得越高的识别率^[24]。

此外，通过计算机仿真时间来衡量成像时间，算法所需时间越短则成像效率越高，算法所需时间越长则可行性越差。

为得到不同信噪比的实测数据，对回波数据添加高斯白噪声。CS重构算法的性能和噪声水平密切相关，而噪声水平是未知的，因此利用能量门限选择噪声单元，然后利用清洁区估计噪声水平^[1]。

阈值的取法，取$\mu =10{{\left\| {{\Theta }^{*}}Y \right\|}_{\infty }}$，平面阈值比原始的线阈值具有更好的压制噪声效果，因为平面阈值考虑到了全局数据，而线阈值仅考虑到当前一列的数据。图 9为不同信噪比条件下OMP,SL0^[25], ExCoV^[26],LB算法对单量测向量ISAR模型的成像结果，其中OMP算法稀疏度冗余设置为40。可以看出，在信噪比较低的情况下，这4种算法都出现了虚假散射点，但SL0算法比OMP算法效果好， ExCoV和LB算法虚假散射点相对较少，由于ExCoV方法仅仅赋予有重要意义的信号元素不同的方差分量，而对剩下的元素赋予单个相同的方差分量^[26]，因此在本文ISAR实测数据散射点系数相差较大条件下，ExCoV方法重构时将强散射点重构出来，而丢失了小散射点，从图中可以明显看出机头部分的散射点丢失了。

图 10为不同信噪比条件下MOMP^[27],MSL0^[28], MLB,MFLB,MILB和MFILB算法对多量测向量ISAR模型的成像结果，其中MOMP算法稀疏度冗余设置为40。可以看出，随着信噪比的降低，MOMP 算法和MSL0算法受信噪比影响严重，低信噪比时出现大量虚假散射点，而本文研究的4种算法受信噪比影响较小，成像效果较好，且优于LB算法，原因在于本文算法使用的平面阈值比LB算法的线阈值具有更好的压制噪声效果，因为平面阈值考虑到了全局数据，而线阈值仅考虑到当前一列的数据。 为进一步定量评价算法的成像性能，下面给出评价指标随信噪比的变化曲线。图 11(a)为不同算法TBR 比较，图 11(b)为不同算法噪声能量比较，图 11(c) 为不同算法相关性比较，图 11(d)为不同算法成像时间比较。通过对比可以看出，相比单量测向量ISAR 模型成像算法以及MOMP算法和MSL0 算法，本文所研究的算法具有更大的TBR，更小的噪声能量， 更大的相关系数，以及更短的成像时间。说明本文算法结果具有更高的识别率，更好的抑制噪声能力，以及较高的可行性。实测数据进一步验证了本文算法的有效性以及成像质量和成像时间的优势。

4 总结

本文针对单量测向量CS-ISAR成像时间长的问题，研究了多量测向量ISAR成像模型，并基于LBI 理论研究了一类多量测向量成像算法，给出了此类算法的整体构架和各种算法的应用条件，然后给出了正则化参数选取方法，在实际中可依据具体的需求选择不同的正则化参数，通过仿真结果表明所提算法具有良好的重构性能及速度优势，基于实测数据的ISAR成像实验也验证了本文算法的有效性，体现出成像时间和成像质量方面的优势。