基于快速SBL的双基地ISAR成像

双基地逆合成孔径雷达(Inverse Synthetic Aperture Radar，ISAR)是发射站和接收站分置的ISAR成像系统，不仅具有双基地雷达的“四抗”特性，而且与单基地ISAR相比，没有几何成像盲区并能获得更多的目标信息[1-3]，日益成为现代雷达技术研究的热点问题之一。

在对目标的观测过程中，通常需要切换雷达波束对多个目标进行多视角观测，容易形成稀疏孔径[4-5]。若直接利用传统的距离-多普勒(Range-Doppler，RD)算法实现双基地ISAR成像，不仅会产生强烈的副瓣和能量泄漏，而且在低信噪比(Signal to Noise Ratio，SNR)条件下存在大量干扰噪声，严重影响成像质量。另外，由于双基地角的存在，其成像分辨率也比相应的单基地雷达低，还容易引起越分辨单元徙动和图像畸变等问题，不利于目标识别[6-7]。

为提高成像质量，国内外不少学者利用ISAR图像的空域稀疏性[8]，将压缩感知(Compressed Sensing，CS)理论[9]应用到了单基地ISAR成像中。文献[10]基于CS理论利用目标像元稀疏先验信息将ISAR图像重构问题转化为基于l1范数稀疏约束最优化问题进行求解，有效改善了成像分辨率和抗噪能力。但在求解时，需要先估算约束参数，且估算误差对成像结果影响较大，有时可能不能获得最优稀疏解。为实现参数自学习，文献[11]提出了基于稀疏贝叶斯学习(Sparse Bayesian Learning，SBL)的ISAR成像算法，利用贝叶斯推理实现参数的迭代自学习，相比于l1范数稀疏约束算法进一步提高了重构质量。但该算法将整个回波矢量化处理，并利用求导的方式求解，涉及大量矩阵求逆等操作，导致数据存储量和运算量较大，运算效率有待提高。另外，考虑到双基地ISAR成像时双基地角的存在[12]，不能直接将单基地成像算法应用到双基地成像中。

基于此，针对稀疏孔径条件下双基地ISAR成像分辨率低、运算时间长等问题，本文提出了一种基于快速SBL的双基地ISAR稀疏孔径高分辨成像算法。首先，对平动补偿后的回波进一步进行多普勒位移补偿，然后将稀疏孔径回波表示为符合CS理论的矩阵形式，再采用分块处理方法利用基于SBL的快速边缘似然函数最大化算法进行图像重构，使得在提高成像质量的同时缩短运算时间，最后将各块回波对应的目标图像合成二维目标图像。通过仿真实验从运算时间和成像性能两个方面验证了算法的有效性和优越性。

1 基于压缩感知的双基地ISAR稀疏孔径成像模型

由于ISAR图像具有很强的空域稀疏性，因此可将CS理论应用到ISAR成像中。本文以文献[6]中的双基地ISAR成像模型为基础，设散射点 P(xP,yP)的散射系数为σP，经过理想的运动补偿后，在tm时刻回波一维距离像可近似表示为

式中，fc为载频中心频率，tp为发射脉冲宽度，μ为调频率，tm为慢时间，

为快时间，θ(tm)和 β(tm)分别为观测时间内的累积转角和双基地角。

1.1 多普勒偏移相位补偿

获得一维距离像后，提取目标的多普勒信息可实现方位成像。根据式(1)，其方位信息包含在相位

中，对相位进行求导，即可得到散射点的多普勒频率：

式中，含有xP的系数项

用于得到目标方位信息，实现方位向高分辨；含有yP的系数项

则可能引起多普勒偏移和图像畸变，由于该项是由相位φ中的第二项引起的，因此需要构造相应的补偿相位将该项补偿掉，构造的补偿相位为

则相位补偿后的一维距离像可表示为

此时，只留下了与散射点多普勒信息有关的相位项，有利于后续方位向成像处理。

1.2 稀疏基矩阵的构造

在应用CS理论时，需要先构造相应的稀疏基矩阵对回波信号进行稀疏表示。假设在距离单元

内共有Q个强散射点，则根据式(4)，此距离单元的回波信号可表示为

式中，aq为第q个散射点的信号复幅度。

假设在观测时间内，全孔径回波信号中共包含L个脉冲视角，其累积转角为Δθ，若将二维成像场景离散化为N个距离单元和M个多普勒单元，则在方位向上散射点q的坐标xq可表示为xq=mΔx(m=1,2,…,M)，其中，

表示每个多普勒单元格在方位向的分辨率，则项

可离散化表示为

式中，l=1,2,…,L。由于收发双站的位置固定以及各脉冲时刻目标到收发双站的距离都是已知的，双基地角可以利用简单的三角几何关系计算出来，且测距精度对双基地角计算的影响几乎可以忽略。因此，将稀疏基矩阵Fall构造为

式中，

当目标平稳运动时，累积转角可近似为均匀变化，则稀疏基矩阵元素可表示为

可以看出，稀疏基矩阵Fall中的元素随着双基地角和累积转角的变化而变化，比标准的傅里叶矩阵更符合回波信号特性，有利于后续的信号重构。

考虑到实际有噪声的存在，将式(5)的回波信号用矩阵形式表示：

式中，Sall表示经预处理后的全孔径二维回波数据，共包含L个回波脉冲，N个距离单元，ε0表示L×N噪声矩阵，Fall表示L×M的稀疏基矩阵，AM×N表示需要求解的ISAR超分辨二维图像，每一个元素对应散射点的复幅度。

1.3 稀疏孔径成像模型

在实际情况下，数据可能存在随机缺失或块缺失两种缺失形式。假设S为融合的有效孔径回波数据，共包含K(K<L)次有效脉冲，则

SK×N=TK×LSallL×N=TFallA+ε=FA+ε

式中，T为有效数据选择矩阵，具有去除缺失孔径，合并有效孔径的作用，F为在稀疏基矩阵Fall中去除缺失孔径对应行后形成的K×M部分稀疏基矩阵，ε为K×N复噪声矩阵，A为需要求解的ISAR超分辨二维图像。为方便求解，将式(9)矢量化表示为

式中，sn，εn和an分别表示第n个距离单元对应的回波数据矢量、噪声矢量和目标图像矢量。

2 快速稀疏贝叶斯学习算法

国内外已有众多学者对基于CS理论的稀疏重构方法进行了研究，并发现与典型的RD算法相比，利用ISAR图像的稀疏特性实现ISAR高分辨成像在提高成像分辨率、降低副瓣和抑制噪声等方面具有明显的优势。但一般的稀疏重构方法仅从优化的角度进行了高分辨重构，有时无法得到最优的稀疏解。SBL算法从统计的角度出发，利用目标的先验信息和贝叶斯推论实现信号重构。由于贝叶斯灵活性强、适用范围广，相比于传统的CS重构算法有一定的优越性[13]，因此本文考虑将SBL算法应用到ISAR成像中，利用ISAR图像的空域稀疏性和先验信息进行建模并实现重构。常用的SBL求解方法主要有直接求导法、期望最大(Expectation Maximization，EM)法和变分贝叶斯(Variational Bayesian，VB)法等，但在求解时涉及大量的矩阵求逆的过程，运算效率较低。在此基础上，Tipping提出了一种快速边缘似然函数最大化法进行求解[14]，并成功应用到回归和分类中，发现该算法不仅在重构精度上与其他算法相差不大，而且大大缩短了运算时间，提高了运算效率。因此，本文将该方法应用到ISAR成像中，提出了一种基于快速SBL的双基地ISAR稀疏孔径高分辨成像算法。

2.1 分块处理模型

若在求解时将全部的二维回波数据矢量化，如式(10)所示，则需要构造复杂的回波数据矢量和对角块稀疏基矩阵，当回波数据较大时，很有可能超过电脑内存空间，且算法执行效率很慢。因此，本文采用拆分的思想，将二维回波矩阵SK×N按距离单元分成小块进行处理，如图1所示。假设每块数据矩阵中包含N1个距离单元的回波数据，则N个距离单元的回波数据一共可分为H=「N/N1⎤块(「 ⎤表示向上取整)，在求解时先逐块数据进行重构，再合成二维图像。

下面只取其中一个数据块回波矩阵进行分析求解。首先根据式(10)将回波矩阵矢量化：

分块处理后，每次重构时处理的数据量从KN×1减小到了KN1×1，相应的对角块矩阵

规模也从KN×MN减小到了KN1×MN1，利于采用SBL算法进行图像重构求解时提高运算效率。由于ISAR成像时雷达回波信号为复数信号，而基于SBL算法的求解过程一般是在实数域推导的，所以需要先把式(11)表示的信号稀疏表示问题由复数域转化为实数域，其转化形式如下：

式中，Re( )和Im( )分别表示实部和虚部。转化之后再利用快速SBL算法求得包含实部和虚部信息的目标矢量

2.2 稀疏先验模型

假设目标图像矢量各像元的实部和虚部ai(i=1,2,…,2MN1)均满足依赖于超参数αi∈α=

的零均值高斯分布，即

则其条件概率密度函数为

假设

是复高斯白噪声，其不同元素εi的虚部和实部分别服从方差为σ2的实高斯分布，有 N(εi|0,σ-2)，则回波信号的条件概率密度函数为

为获得高斯先验的共轭特性，对超参数αi和σ-2施加同高斯对偶的Gamma先验分布：

p(σ-2)=Gamma(σ-2|c,d)

式中，Gamma(α|a,b)=Γ(a)-1baαa-1e-bα， Γ(a)=

tα-1e-1d t。

可以看出，目标图像矢量

实际上是通过超参数α实现控制的，而超参数α又由系数a,b约束，所以对目标图像的稀疏促进作用可以看作是一个两层稀疏先验模型，如图2所示。

2.3 快速边缘似然函数最大化求解

为缩短运算时间，采取文献[14]中的快速边缘似然函数最大化方法进行求解。在贝叶斯推理中结合目标先验信息和似然函数，得到目标矢量

的后验概率为[13]

可以看出目标矢量

的后验概率也服从高斯分布，即

且

式中，Λ=diag(α1,α2,…,α2MN1)表示对角元素由α1,α2,…,α2MN1组成的对角矩阵。此时的均值μ即为对应的目标矢量估计值

可通过最大化

得到α,σ2的估计值[13]。由于

则使

的最优解同时也是使

取得最大值的解，为方便计算，在对数域构造代价函数：

式中，

当只考虑

中第i个基原子Fi时，其对应系数的超参数为αi，为使代价函数最大化，将矩阵C进行分解，可得到

式中，C-i表示C中与基向量Fi无关的部分。根据矩阵行列式及求逆计算规律展开式(20)，可得

式中，α-i表示从α中除去第i个超参数αi后剩下的超参数集合，且有

可以看出，此时代价函数可以分解为两部分，第一部分L(α-i)为

中除去第i列基Fi的其他列基对代价函数的影响，第二部分l(αi)为只考虑基Fi对代价函数的影响。为得到代价函数最大值，先用L(α)对αi求一阶偏导数，可得

对该式可分为两种情况进行分析：

1) 当

时，式(25)有唯一解

为进一步分析函数在该驻点的性质，用L(α)对αi求二阶偏导数，可得

式(27)表示，代价函数L(α)在唯一解

处取得极大值。

2) 当

时，式(25)结果总大于零，说明L(α)是递增函数，则当αi=+∞时取得极大值。

由于直接计算qi,gi比较复杂，故考虑引入两个中间变量Gi和Qi以简化运算，有

式中，B=σ-2I，则

当αi=+∞时，有gi=Gi，qi=Qi，因此通过计算Gi和Qi即可得到

中基原子Fi对应的gi和qi，从而避免了大量的复杂计算。

2.4 算法流程图

由上述分析可知，对基矩阵

中的第i列基Fi，当满足

且αi<+∞时，对其超参数αi进行更新，有

当满足

且αi=+∞时，将Fi添加到被选中的原子中；当满足

时，则将Fi从被选中的原子中删除，并令αi=+∞。通过该原则按顺序对基函数进行添加、更新或删除并对其对应的超参数进行更新，就能实现快速SBL算法求解，求解过程如图3所示，具体的实现步骤可描述为：

1) 输入回波矢量

和稀疏基对角块矩阵

初始化迭代次数n=1；

2) 将矩阵

中的每一列基作为一个原子，将每个原子与回波矢量

作内积，选择内积值最大的原子作为被选中的第一个原子Fi，利用式

和

进行初始化，其中

表示

的标准差；

3) 利用式(18)和式(19)计算∑,μ，再根据式(28)和式(29)计算出

中每个原子对应的Gi和Qi；

4) 利用式(30)和式(31)计算出每个原子对应的gi和qi以及

值，分3种情况进行处理：

① 若θi>0且αi=+∞，则将Fi添加到被选中的原子中；

② 若θi>0且αi<+∞，则对超参数αi进行更新，有

③ 若θi<0，则将Fi从被选中的原子中删除，并设置αi=+∞；

5)利用

更新σ2[13]，得到新的

和σ2后，再利用式(19)和式(18)更新∑,μ，利用式(28)和式(29)计算新的Gi和Qi；

6) 当迭代次数达到设定值或估计值满

时终止迭代，此时的μ值即为目标矢量的估计值

否则置n=n+1并返回第2)步继续执行；

7) 将所求目标矢量估计值

中的实部与相应的虚部相结合变为复数矢量信号，再将其转化为二维矩阵形式，得到二维目标图像

综上所述，基于快速SBL算法的双基地ISAR稀疏孔径高分辨成像过程可描述为：首先，为解决双基地角时变带来的多普勒频移问题，构造补偿相位，得到补偿后的二维回波数据Sall，然后利用数据选择矩阵T，得到方位向稀疏孔径二维回波数据SK×N，再按距离单元分为H块，每块只包含N1个距离单元的回波数据，将回波表示为

并将每块回波数据进行矢量化并转化到实数域，然后利用快速SBL算法分别对每块回波的目标图像矢量进行求解，得到每块的目标图像估计值，最后合成整个二维图像，得到稀疏孔径目标二维图像估计值

3 仿真实验与分析

本节主要从运算时间和成像性能这两个方面对算法性能进行仿真验证。为验证本文算法在缩短运算时间方面的优越性，利用快速边缘似然函数最大化法以及直接求导法这两种SBL求解方法对一维稀疏信号进行重构，通过比较重构的均方根误差和重构时间，体现算法在运算时间方面的优越性。为验证算法在成像性能方面的优越性，利用典型的RD算法、文献[10]中提出的基于l1范数稀疏约束最优化算法以及本文算法对双基地ISAR稀疏孔径回波进行目标成像，通过对比重构图像的目标背景比(Target-to-Background Ratio，TBR)[15]和图像熵，体现算法在成像性能方面的优越性。

3.1 一维稀疏信号重构性能仿真

为验证算法在运行时间方面的优越性，选取简单的线性调频信号作为原始信号进行重构。假设信号为x(t)=cos(2π(1 000t+(10t2)/2))，采样频率为400 Hz，信号长度为N=500，对原始信号进行随机采样，稀疏基为傅里叶基，利用本文提到的快速边缘似然函数最大化算法和直接求导法对SBL模型进行求解，总迭代次数均设为1 000次，收敛条件均设为

在观测点数为M=250时利用快速算法的重构结果如图4所示。改变观测点数，采用运算时间和均方根误差(RMSE=

衡量信号的重构性能，如表1所示。从表中数据可以看出，在不同的观测点数情况下，两种求解方法的重构精度基本一致，但在重构时间上有较大差别，与直接求导方法相比，快速边缘似然函数最大化方法在保证重构精度的同时大大提高了运算速度，体现了其在运算时间方面的优越性。

3.2 双基地ISAR稀疏孔径二维成像仿真

为验证算法的成像性能，利用典型的RD算法、基于l1范数稀疏约束最优化算法及本文算法对双基地ISAR稀疏孔径回波进行目标图像重构。双基地ISAR仿真场景如图5所示，假设双基地基线长度为400 km，目标在300 km的高度以3 km/s的速度匀速运动，成像起点为水平距离距接收站雷达靠右70 km处。目标的散射点模型如图6所示，该模型由96个散射点组成，成像仿真参数设置如表2所示，累积脉冲数为400，截取其中400个距离单元数据进行处理，在不同的SNR和孔径缺失条件下，利用3种算法进行成像，其中本文算法以每20个距离单元的回波数据为一个独立的小块进行处理。

3.2.1 不同孔径缺失情形下成像性能验证

为验证算法在不同孔径缺失情形下的成像性能，在SNR为5 dB条件下，通过改变孔径缺失情形，将本文算法与基于l1范数稀疏约束最优化算法和RD算法的成像结果进行对比。4种孔径缺失情形分别是50%数据随机缺失、75%数据随机缺失、50%数据块缺失和75%数据块缺失，利用3种算法进行成像的结果如图7所示，其中，从左至右每列分别对应RD算法、基于l1范数稀疏约束最优化算法和本文算法的成像结果。

从成像结果可以看出，在孔径稀疏的情况下，采用RD算法直接成像会有大量的能量泄漏并产生散焦现象。相比之下，采用基于l1范数稀疏约束最优化算法成像能有效提高成像质量，能量泄漏问题和图像散焦问题得到改善，但仍存在少量噪声没有得到抑制。另外，当孔径缺失数量较多时(如75%数据块缺失的情况)，利用该方法重构图像只能将目标的基本轮廓重构出来，而无法有效地区分各个散射点，成像质量迅速下降。采用本文算法进行高分辨成像时，不仅能有效抑制噪声，而且在数据缺失较多时仍能高质量地重构出目标图像。从表3中两种稀疏重构方法的评价指标也能看出，在同一孔径缺失情形下，相比于基于l1范数稀疏约束最优化算法，本文算法重构图像的TBR值大、图像熵值小，且运算时间短，说明了所提算法能在提高成像质量的同时减少运算时间，体现了算法的优越性。

3.2.2 不同SNR条件下成像性能验证

为验证算法在不同SNR条件下的成像性能，在方位向孔径随机缺失比为50%的情况下，利用基于l1范数稀疏约束最优化算法和本文算法在不同SNR条件下进行稀疏孔径成像，并画出目标图像的TBR和熵值随SNR变化的曲线，如图8所示。从变化曲线图可以看出，在同一SNR条件下，相较于基于l1范数稀疏约束最优化算法，利用本文算法所得图像的TBR值大、熵值小，说明本文算法能实现更高质量的图像重构。另外，随着SNR降低，虽然两种算法的成像结果的TBR值都有所减小、熵值有所增大，但本文算法变化较小，说明SNR对算法的影响较小，在低SNR条件下本文算法仍然适用。

4 结束语

本文提出的基于快速SBL的双基地ISAR稀疏孔径高分辨成像算法通过快速边缘似然函数最大化方法进行SBL问题求解，相比于直接求导方法在保证重构精度的同时大大缩短了运算时间。另外，为减少每次数据处理量和存储量，采用分块处理的思想，将整个回波矩阵分为若干个小块进行处理，进一步提高了运算效率。仿真实验表明，在孔径缺失数据较多以及低SNR条件下，本文算法仍能得到高质量目标图像。

[1] 保铮,邢孟道,王彤. 雷达成像技术[M]. 北京: 电子工业出版社, 2005.

[2] 杨振起,张永顺,骆永军. 双(多)基地雷达系统[M]. 北京: 国防工业出版社, 1998.

[3] CATALDO D, MARTORELLA M. Super-Resolution for Bistatic Distortion Mitigation[C]∥ IEEE Radar Conference, Philadelphia, PA: IEEE, 2016:1-6.

[4] 李少东,陈文峰,杨军,等. 稀疏孔径下的运动补偿及快速超分辨成像方法[J]. 电子学报, 2017,45(2):291-299.

[5] ZHANG Shunsheng, ZHANG Wei, ZONG Zhulin, et al. High-Resolution Bistatic ISAR Imaging Based on Two-Dimensional Compressed Sensing[J]. IEEE Trans on Antennas and Propagation, 2015, 63(5):2098-2111.

[6] 郭宝锋,尚朝轩,王俊岭,等. 双基地角时变下的逆合成孔径雷达越分辨单元徙动校正算法[J]. 物理学报, 2014, 63(23):238406.

[7] GUO Baofeng, WANG Junling, GAO Meiguo, et al. Research on Spatial-Variant Property of Bistatic ISAR Imaging Plane of Space Target[J]. Chinese Physics B, 2015, 24(4):048402.

[8] XU Gang, XING Mengdao, XIA Xianggen, et al. High-Resolution Inverse Synthetic Aperture Radar Imaging and Scaling with Sparse Aperture[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015, 8(8):4010-4027.

[9] 钱宇雷. 基于压缩感知的雷达成像方法研究[D]. 南京: 南京航空航天大学, 2016.

[10] ZHANG Lei, QIAO Zhijun, XING Mengdao, et al. High-Resolution ISAR Imaging by Exploiting Sparse Apertures[J]. IEEE Trans on Antennas and Propagation, 2012, 60(2):997-1008.

[11] LIU Hongchao, JIU Bo, LIU Hongwei, et al. Superresolution ISAR Imaging Based on Sparse Bayesian Learning[J]. IEEE Trans on Geoscience and Remote Sensing, 2014, 52(8):5005-5013.

[12] BAE J, KAE B, LEE S, et al. Bistatic ISAR Image Reconstruction Using Sparse-Recovery Interpolation of Missing Data[J]. IEEE Trans on Aerospace and Electronic Systems, 2016, 52(3):1155-1167.

[13] TIPPING M E. Sparse Bayesian Learning and the Relevance Vector Machine[J]. Journal of Machine Learning Research, 2001,1:211-244.

[14] TIPPING M E, FAUL A C. Fast Marginal Likelihood Maximisation for Sparse Bayesian Models[C]∥Ninth International Workshop on Artificial Intelligence and Statistics, Key West, FL:[s.n.], 2003:1-13.

[15] ZHANG Lei, WANG Hongxian, QIAO Zhijun. Resolution Enhancement for ISAR Imaging via Improved Statistical Compressive Sensing[J]. EURASIP Journal on Advances in Signal Processing, 2016, 2016[80]:1-19.