1.   引言

由于具有全天时、全天候、远作用距离、高分辨率等独特优势，逆合成孔径雷达(Inverse Synthetic Aperture Radar, ISAR)广泛用于高分辨成像、目标特征提取与分类识别等领域^[1,2]。通过发射稀疏步进调频信号^[3]，ISAR可以在获得高距离分辨率的同时减小接收机瞬时带宽，从而降低硬件要求。同时，通过设计子脉冲载频序列可以缩短观测时间并提高抗干扰能力。但是，稀疏步进调频信号对目标径向运动十分敏感，运动补偿后的残余径向运动将导致ISAR像距离和方位维的严重散焦。同时，非均匀分布的子脉冲序列使传统基于傅里叶变换的距离像合成方法失效。此外，由于子脉冲能量有限，在远距离、小目标观测时回波信噪比较低，实际成像中，一般认为0 dB回波信噪比为低信噪比^[4]。因此，如何在低信噪比等复杂观测环境下，实现稀疏步进调频信号的运动补偿及聚焦成像近年来受到了雷达成像领域的广泛关注^[5,6]。

稀疏步进调频信号的平动包含：(1)子脉冲间平动；(2)脉冲串间平动。针对步进调频信号，主要采用最小平均距离像熵^[7]、数字拉伸^[8]等运动补偿方法。针对稀疏步进调频信号，文献[9]基于补充编码对消方法实现对目标速度的精确估计。文献[10]设计了联合熵代价函数，并采用粒子群滤波(Particle Swarm Optimization, PSO)方法估计子脉冲及脉冲串间的平动。但是，上述方法存在过程复杂、精度不高等问题。

在高分辨1维距离像(High-Resolution Range Profile, HRRP)合成方面，由于稀疏步进调频信号在解线频调处理后往往对应非均匀的频率网格，因此无法采用傅里叶变换成像。针对该问题，基于稀疏信号重构^[11,12]的HRRP合成方法受到了广泛关注。该类方法主要包括数值优化和稀疏贝叶斯学习方法。其中，数值优化的典型方法有正交匹配追踪(Orthogonal Matching Pursuit, OMP)算法^[10,13]、L1范数优化^[14]等。但是，上述方法在低信噪比条件下重构误差较大。基于稀疏贝叶斯学习的方法首先对散射点引入稀疏先验并构建概率图模型，进而采用近似推断、采样等方法进行求解。文献[15]令散射点服从Laplace分布，进而采用梯度下降法对随机相位及HRRP进行求解。但是，Laplace分布仅有一个尺度参数，因此较难单独调整散射点向量中的某些元素以实现稀疏重构。

本文通过构造参数化字典，将稀疏步进调频信号的低信噪比ISAR成像问题转换为运动参数估计与ISAR像重构联合估计问题。同时针对现有运动参数估计方法在低信噪比时误差大，容易陷入局部最优解等不足，提出遗传算法与基于Gamma-Gauss先验的变分贝叶斯推断(Variational Bayesian Inference, VBI)相结合的稀疏步进调频信号低信噪比ISAR成像方法。该方法首先根据搜索区间内的运动参数构造字典，进而采用VBI进行低信噪比HRRP合成，最终以所得2维ISAR像的图像熵最小为准则对运动参数种群进行更新。本方法能够有效解决低信噪比下HRRP重构误差对运动参数估计的不利影响，同时能够通过种群更新跳出局部最优值，最终获得目标运动参数的精确估计与高分辨聚焦成像。仿真及实测数据的处理结果验证了所提算法的有效性。

2.   信号模型

2.1   回波模型

通过设定载频序列，可以从完整的步进调频信号中抽取得到稀疏步进调频信号。令完整的步进调频信号包含$N$个脉冲，并从中依次抽取$M$个脉冲($M < N$)。设稀疏步进调频信号的载频序列为${f_{sm}} =$${f_{\rm{c}}} + g\left( m \right)\Delta f$，$\Delta f$对应脉冲带宽，$g$为$\left[ {0:N - 1} \right]$的一个子集，则稀疏步进调频信号的带宽为$B = N\Delta f$。分析可知，其相干积累时间小于完整的步进调频信号。同时，可根据环境信息设定载频序列，从而具备一定的抗干扰能力。

假设雷达共发射$K$组稀疏步进调频信号，则第$k$($k = 1,2, \cdots ,K$)组稀疏步进调频信号为

其中，$t = \hat t{\rm{ + }}m{T_{\rm{R}}} + kM{T_{\rm{R}}}$($m = 1,2, \cdots ,M$)为全时间，$\hat t$为快时间，${\rm{rect}}\left( u \right)$表示矩形窗，当$\left| u \right| \le {1 / 2}$时，${\rm{rect}}\left( u \right)$为1，否则${\rm{rect}}\left( u \right)$为0, $\gamma$为调频率，${T_{\rm{p}}}$和${T_{\rm{R}}}$分别表示脉冲宽度和脉冲重复周期。

假设目标包含$P$个散射点，第$p$($p = 1,2, \cdots ,P$)个散射点的后向散射系数为${\sigma _p}$。基于Stop-Go模型^[16]，令一个脉冲内散射点$p$的时延${\tau _p}\left( t \right)$不变，即有${\tau _p}\left( t \right) \approx {\tau _p}\left( {{t_{m,k}}} \right)$, ${t_{m,k}} = m{T_{\rm{R}}} + kN{T_{\rm{R}}}$。则对于该散射点，第$k$组稀疏步进调频信号的第$m$个子脉冲回波为

其中，${\tau _p}\left( {{t_{m,k}}} \right) = {{2{R_p}\left( {{t_{m,k}}} \right)} / {\rm{c}}}$, ${R_p}\left( {{t_{m,k}}} \right)$为第$p$个散射点到雷达的瞬时斜距，${\rm{c}}$为光速，$\varepsilon \left( {\hat t} \right)$为加性噪声。当成像积累角${\theta _{m,k}}$较小时，${R_p}\left( {{t_{m,k}}} \right) = R\left( {{t_{m,k}}} \right) +$${x_p}\sin {\theta _{m,k}} + {y_p}$，其中$R\left( {{t_{m,k}}} \right)$为雷达到参考点的瞬时斜距，$\left( {{x_p},{y_p}} \right)$为散射点$p$在成像平面上的坐标。对于高速运动目标，${t_{m,k}}$时刻目标参考点到雷达的距离$R\left( {{t_{m,k}}} \right) = {r_{\rm{R}}} + {v_{\rm{R}}}{t_{m,k}} + {1 / 2}{a_{\rm{R}}}t_{m,k}^2$，其中${r_{\rm{R}}}$为初始时刻目标参考点到雷达的距离，${v_{\rm{R}}}$为目标径向速度，${a_{\rm{R}}}$为目标径向加速度。参考信号时延为${\tau _{{\rm{ref}}}}\left( {{t_{m,k}}} \right) = {{2{R_{{\rm{ref}}}}\left( {{t_{m,k}}} \right)} / {\rm{c}}}$, ${R_{{\rm{ref}}}}\left( {{t_{m,k}}} \right) = {r_{\rm{R}}} + {\hat v_{\rm{R}}}{t_{m,k}} +$${1 / 2}{\hat a_{\rm{R}}}t_{m,k}^2$。通常，${v_{\rm{R}}}$, ${a_{\rm{R}}}$的粗估计${\hat v_{\rm{R}}}$, ${\hat a_{\rm{R}}}$可在目标跟踪阶段获得。对回波进行解线频调处理并令$\hat f = \gamma \left( {\hat t - {{2R\left( {{t_{m,k}}} \right)} / {\rm{c}}}} \right)$，则信号化简为

其中，$\Delta R = {x_p}\sin {\theta _{m,k}} + {y_p}$, ${\varPhi _{\rm{P}}}$, ${\varPhi _{\rm{B}}}$分别为子脉冲间平动引起的相位误差、脉冲串间平动引起的相位误差。${\varPhi _{\rm{P}}}$, ${\varPhi _{\rm{B}}}$表达式为

其中，$\Delta {v_{\rm{R}}}$为剩余速度，$\Delta {v_{\rm{R}}} = {v_{\rm{R}}} - {\hat v_{\rm{R}}}$, $\Delta {a_{\rm{R}}}$为剩余加速度，$\Delta {a_{\rm{R}}} = {a_{\rm{R}}} - {\hat a_{\rm{R}}}$, ${\varPhi _1}$为$m$的2次相位项，距离像合成时会造成主瓣展宽；${\varPhi _2}$为$m$的3次相位项，合成距离像时会造成非对称旁瓣，该项数量级一般为10^–4，因此可忽略；${\varPhi _3}$包含$m$, $k$的耦合项，会造成距离像偏移，造成包络弯曲；${\varPhi _4}$为$k$的1次相位项，会造成图像的方位向偏移，可忽略；${\varPhi _5}$是$k$的2次相位项，方位向脉压时会造成方位像主瓣展宽。3.3节重点分析${\varPhi _1}$, ${\varPhi _3}$, ${\varPhi _5}$对成像的影响。

为构建基于稀疏表示的HRRP重构算法，将回波改写为离散形式。假设每个脉冲采样点数为${N_{\rm{r}}}$，并令$\bar L = M \cdot {N_{\rm{r}}}$，则第$k$组回波可表示为${{\boldsymbol{s}}_k} =$$\left[ {{{\boldsymbol{s}}_{0,k}},\cdots,{{\boldsymbol{s}}_{1,k}}, \cdots ,{{\boldsymbol{s}}_{M - 1,}}_k} \right]_{_{1 \times \bar L}}^{\rm{T}}$，其中，${{\boldsymbol{s}}_{m,k}} = \left[ {s_{0,m,k,}}{\rm{ }}\cdots\right.$$\left.{s_{1,m,k,}} \cdots{s_{{N_{\rm{r}}} - 1,m,k}} \right]_{1 \times {N_{\rm{r}}}}^{}$, ${s_{{n_{\rm{r}}},m,k}} = \exp \left( {\rm{j}}\dfrac{{{\rm{4\pi }}}}{{\rm{c}}}\left( {{f_{sm}} + \dfrac{{{n_r}}}{{{N_{\rm{r}}}}}\Delta f} \right)\right.$$\left.\left( \Delta R +\Delta {v_{\rm{R}}}{t_{m,k}} + \dfrac{1}{2}\Delta {a_{\rm{R}}}t_{m,k}^2 \right) \right) + \varepsilon \left( {{n_{\rm{r}}}} \right)$。

将目标剩余径向运动参数$\left( {\Delta {v_{\rm{R}}},\Delta {a_{\rm{R}}}} \right)$引入字典，则稀疏观测模型为

其中，${{\boldsymbol{D}}_k}\left( {\Delta {v_{\rm{R}}},\Delta {a_{\rm{R}}}} \right) \in {\mathbb{C}^{\bar L \times L}}$为第$k$次回波对应的字典矩阵，$L = N \cdot {N_{\rm{r}}}$, ${{\boldsymbol{\theta}} _k} \in {\mathbb{C}^{L \times 1}}$为第$k$次回波对应的HRRP，${\boldsymbol{n}}$为噪声向量。定义${\boldsymbol{d}}_l^k$为字典矩阵${{\boldsymbol{D}}_k}\left( {\Delta {v_{\rm{R}}},\Delta {a_{\rm{R}}}} \right)$的第 $l$列，则 ${\boldsymbol{d}}_l^k$表达式为

其中，$\odot$表示内积，${{\boldsymbol{f}}_l}$为字典${\boldsymbol{F}}$第$l$列，${\boldsymbol{F}} =$$\left[ {{{\boldsymbol{F}}_0}, {{\boldsymbol{F}}_1},\cdots ,{{\boldsymbol{F}}_{M - 1}}} \right]_{\bar L \times L}^{\rm{T}}$, ${{\boldsymbol{g}}_k} = \left[ {{{\boldsymbol{g}}_{0,k}}, {{\boldsymbol{g}}_{1,k}},\cdots ,{{\boldsymbol{g}}_{M - 1,k}}} \right]_{\bar L \times 1}^{\rm{T}}$, ${{\boldsymbol{h}}_k} = \left[ {{{\boldsymbol{h}}_{0,k}},{{\boldsymbol{h}}_{1,k}}, \cdots ,{{\boldsymbol{h}}_{M - 1,k}}} \right]_{\bar L \times 1}^{\rm{T}}$。

对于第$k$组信号第$m$个子脉冲来说

2.2   贝叶斯模型

由于基于稀疏贝叶斯学习的信号重构方法通常对实数矩阵进行操作，因此将式(5)改写为

其中，${\boldsymbol{y}} = {\left[ {{\rm{Re}} \left( {{{\boldsymbol{s}}_k}} \right),{\rm{Im}} \left( {{{\boldsymbol{s}}_k}} \right)} \right]^{\rm{T}}}$, ${\boldsymbol{\theta}} = {\left[ {{\rm{Re}} \left( {{{\boldsymbol{\theta}} _k}} \right),{\rm{Im}} \left( {{{\boldsymbol{\theta}} _k}} \right)} \right]^{\rm{T}}}$, ${\boldsymbol{\varPhi}} = \left[ {\begin{array}{*{20}{c}} {{\rm{Re}} \left( {{{\boldsymbol{D}}_k}} \right)}&{ - {\rm{Im}} \left( {{{\boldsymbol{D}}_k}} \right)} \\ {{\rm{Im}} \left( {{{\boldsymbol{D}}_k}} \right)}&{{\rm{Re}} \left( {{{\boldsymbol{D}}_k}} \right)} \end{array}} \right]$, ${\boldsymbol{\varepsilon}} = {\left[ {{\rm{Re}} \left( {\boldsymbol{n}} \right),{\rm{Im}} \left( {\boldsymbol{n}} \right)} \right]^{\rm{T}}}$。

针对模型式(10)，通过引入Gamma-Gauss先验构建概率图模型。假设噪声${\boldsymbol{\varepsilon}}$为0均值的高斯白噪声，此时回波${\boldsymbol{y}}$的概率分布同样为高斯分布，${\boldsymbol{\varepsilon}}$与${\boldsymbol{y}}$概率分布为

其中，$\alpha$为噪声精度。对稀疏向量${\boldsymbol{\theta}}$引入Gamma-Gauss先验

其中，精度矩阵${\boldsymbol{\varLambda}} = {\rm{diag}}\left( {{\lambda _1},{\lambda _2},\cdots,{\lambda _D}} \right)$为$D$维对角阵，$D = 2L$。精度参数${\lambda _d}$($d = 1,2,\cdots,D$), $\alpha$服从Gamma分布

对应的概率图模型如图1所示。根据概率图模型，联合分布可表示为所有分布的乘积，即

而随机变量的后验分布可表示为各变量的联合分布除以边缘分布$p\left( {\boldsymbol{y}} \right)$，从而得到

直接计算模型后验分布十分困难，因此利用VBI算法求解近似后验分布，从而实现低信噪比条件下的HRRP合成。

图 1  概率图模型

下载: 全尺寸图片 幻灯片

3.   基于VBI和遗传算法的运动目标ISAR成像算法

3.1   目标运动参数估计及高分辨成像方法

机动目标稀疏步进调频信号成像的核心问题是低信噪比下的有效的HRRP合成及精确的剩余运动补偿。根据式(10)模型，可通过求解如下目标函数获得目标散射点分布与剩余运动参数的准确估计

其中，$E(\cdot)$为图像熵，$I\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$为补偿径向运动参数$\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$后的2维ISAR图像。由于式(16)代价函数包含较多的局部最优值，采用梯度下降法难以保证全局最优解，因此采用遗传算法求解${P_1}$^[17]问题。算法主要步骤如下：

步骤1　初始化。令${\rm{iter1}} = 1$，随机初始化$\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$种群并设置终止次数${{{G}}_{\rm{1}}}$；

步骤2　对每组$\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$构造相应的${\boldsymbol{D}}_k^{\left( {{\rm{iter1}}} \right)}\cdot$$\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$，并重写信号模型为

根据所构建的概率图模型求解${\boldsymbol{\theta}} _k^{\left( {{\rm{iter1}}} \right)}$，并构建距离像矩阵${\boldsymbol{s}}_{\rm{r}}^{\left( {{\rm{iter1}}} \right)}\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right) = \left[ {\boldsymbol{\theta}} _1^{\left( {{\rm{iter1}}} \right)},{\boldsymbol{\theta}} _2^{\left( {{\rm{iter1}}} \right)}, \cdots ,$${\boldsymbol{\theta}} _K^{\left( {{\rm{iter1}}} \right)} \right]$，继而获得相应的ISAR 2维图像${I^{\left( {{\rm{iter1}}} \right)}}$$\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$。计算图像熵

其中，${I'_{{n_{\rm{r}}},{n_{\rm{a}}}}}\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right) = {{I_{{n_{\rm{r}}},{n_{\rm{a}}}}}\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)} \Big/$${\displaystyle\sum\nolimits_{{n_{\rm{a}}}}^{{N_{\rm{a}}}} {\sum\nolimits_{{n_{\rm{r}}}}^{{N_{\rm{r}}}} {{I_{{n_{\rm{r}}},{n_{\rm{a}}}}}\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)} } }$为归一化图像；

步骤3　保留图像熵较小的$\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$个体并更新$\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$种群，令${\rm{iter1 = iter1 + 1}}$；

步骤4　循环步骤2和步骤3直至达到循环次数${G_{\rm{1}}}$。

根据上述流程，步骤2的关键是准确重构HRRP，重构算法将在3.2节详细讨论。步骤3将ISAR图像熵作为遗传算法适应度值，通过基因选择、基因交叉及基因变异3个步骤更新种群。基因交叉及基因变异可产生新种群，增加了种群的多样性并帮助种群跳出局部最优值。

3.2   基于VBI算法的高分辨距离像重构

由于本文采用共轭先验进行概率建模，因此可以直接采用VBI算法^[18]求得后验分布的闭式解。假设${\boldsymbol{X}}$表示观测数据，${\boldsymbol{w}}$表示随机变量的集合，则${w_j}$的近似后验分布的表达式为

其中，${\rm{const}}$为保证后验分布归一化的常数。接下来，根据式(19)以及2.2节构建的概率模型，利用VBI算法分别求解变量${\boldsymbol{\theta}}$, $\alpha$, ${\boldsymbol{\varLambda}}$的近似后验分布以实现稀疏重构。根据平均场(mean field)假设，令式(15)为${\boldsymbol{\theta}}$, $\alpha$, ${\boldsymbol{\varLambda}}$后验相乘的形式

下面给出${\boldsymbol{\theta}}$求解步骤：

步骤1　初始化，给定${\boldsymbol{\varPhi}}$，初始化${v_1}$, ${v_2}$, ${v_3}$, ${v_4}$, ${\boldsymbol{\varLambda}}$, $\alpha$，终止阈值${\eta _1}$，终止次数${{{G}}_{\rm{2}}}$；

步骤2　更新变量$\alpha$，根据式(19)，其近似后验分布为Gamma分布

其中，$a' \,=\, {v_3} \,+\, {{\bar L} / 2}$, $b' \,=\, {v_4} \,+\, \dfrac{1}{2}{E_{q\left( {\boldsymbol{\theta}} \right)}}$ $\cdot\left[ {{{\left( {{\boldsymbol{y}} - {\boldsymbol{\varPhi}} {\boldsymbol{\theta}} } \right)}^{\rm{T}}}\left( {{\boldsymbol{y}} - {\boldsymbol{\varPhi}} {\boldsymbol{\theta}} } \right)} \right]$，$\bar L$为${\boldsymbol{y}}$的长度。此时，$\alpha$的期望为

步骤3　更新变量${\boldsymbol{\varLambda}}$，其第$d$个元素${\lambda _d}$的近似后验分布为Gamma分布

其中，${e'_d} = {v_1} + \dfrac{1}{2}$, ${f'_d} = {v_2} + \dfrac{1}{2}{{\rm{E}}_{q( {\boldsymbol{\theta}} )}}\left[ {{\boldsymbol{\theta}} _d^2} \right]$。${\lambda _d}$的期望为

步骤4　更新变量${\boldsymbol{\theta}}$，其近似后验分布为高斯分布

其中，${\boldsymbol{\varSigma}} ' = {\left( {{{\rm{E}}_{q\left( {\boldsymbol{\varLambda}} \right)}}\left[ {\boldsymbol{\varLambda}} \right] + {{\rm{E}}_{q\left( \alpha \right)}}\left[ \alpha \right]{{\boldsymbol{\varPhi}} ^{\rm{T}}}{\boldsymbol{\varPhi}} } \right)^{ - 1}}$, ${\boldsymbol{\mu}} ' =$${{\rm{E}}_{q\left( \alpha \right)}}\left[ \alpha \right]{\boldsymbol{\varSigma}} '{{\boldsymbol{\varPhi}} ^{\rm{T}}}{\boldsymbol{y}}$。此时

步骤5　重复步骤步骤2—步骤4，当$\hat {\boldsymbol{\theta}}$相对于前一次估计的相对变化量不超过阈值${\eta _1}$时，停止迭代。

综上所述，基于VBI和遗传算法的稀疏步进调频信号高分辨ISAR成像算法的伪代码如表1所示。

表 1  所提算法伪代码

算法：基于稀疏步进调频信号的低信噪比ISAR成像
(1) 初始化种群$\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$, ${v_1}$, ${v_2}$, ${v_3}$, ${v_4}$, ${\boldsymbol{\varLambda}}$, $\alpha$, ${{{G}}_{\rm{1}}}$, ${{{G}}_{\rm{2}}}$, ${\eta _1}$。
(2) For iter1=1:${{{G}}_{\rm{1}}}$
(a) For k=1:${{K}}$
构造${\boldsymbol{D}}_k^{\left( {{\rm{iter1}}} \right)}\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$；
For iter2=1:${{{G}}_{\rm{2}}}$
利用式(22)更新$\alpha$；
利用式(24)更新${\boldsymbol{\varLambda}}$；
利用式(26)更新${\boldsymbol{\theta}}$；
若$\hat {\boldsymbol{\theta}} _k^{\left( {{\rm{iter1}}} \right)}$相对于前一次估计的变化量小于${\eta _1}$则停止循环；
End
End
(b) 构建距离像矩阵${\boldsymbol{s} }_{\rm{r} }^{\left( { {\rm{iter1} } } \right)}\left( {\Delta { {\hat v}_{\rm{R} } },\Delta { {\hat a}_{\rm{R} } } } \right) = \left[ { {\boldsymbol{\theta} } _{\rm{1} }^{\left( { {\rm{iter1} } } \right)},{\boldsymbol{\theta} } _{\rm{2} }^{\left( { {\rm{iter1} } } \right)}, \cdots,{\boldsymbol{\theta} } _K^{\left( { {\rm{iter1} } } \right)} } \right]$，得${I^{\left( {{\rm{iter1}}} \right)}}\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$；
(c) 根据式(18)计算图像熵；
(d) 保留图像熵较小的个体并更新种群；
(e) 判断是否达到循环次数${{{G}}_{\rm{1}}}$；
End
(3) 实现高分辨ISAR成像。

下载: 导出CSV 

| 显示表格

3.3   算法分析

下面将根据式(4)分析对$\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$的精度要求。对于HRRP合成来说，2次相位项${{{\varPhi}} _1}$将造成距离像展宽，因此要求脉冲串相干积累时间内${{{\varPhi}} _1}$变化小于${{\rm{\pi }} / {\rm{2}}}$，则

在实际中，式(27b)通常可以满足，因此忽略$\Delta \hat a$对HRRP合成的影响。为了消除1阶相位项${{{\varPhi}} _3}$造成的包络弯曲，要求成像观测时间内剩余速度和剩余加速度引起的包络偏移不超过半个距离分辨单元^[19]，则

值得注意的是，$\Delta \hat v$在脉冲串间产生线性相位引起距离像偏移与弯曲，而$\Delta \hat a$导致方位主瓣展宽与散焦。

为获得良好的方位聚焦效果，要求成像观测时间内剩余加速度引起的方位像峰值降低不超过3 dB，则

因此，方位聚焦需要更高的加速度估计精度。最终，为获得聚焦效果良好的图像，目标剩余速度和剩余加速度应满足式(28a)及式(29)。后续实验表明，所提算法可满足该估计精度要求。

4.   仿真实验与分析

4.1   蒙特卡罗实验

为定性验证VBI算法的优越性，将其与典型数值优化、稀疏贝叶斯学习方法进行比较，利用相变图比较OMP算法^[10]、梯度下降(Gradient Descent, GD)算法^[15]、VBI算法的性能。相变图通过多次蒙特卡罗实验进行信号重构，采用归一化均方误差(Normalized Mean Square Error, NMSE)分析重构结果准确性^[20]。相变图的相变区域分为不可重建和可重建区域，能够对算法重构性能进行直观描述。首先构建SNR和采样比组成的2维相变平面，假设稀疏向量${\boldsymbol{\theta}}$长度为1024，如图2所示，非0元素个数为18，幅值为–1或1。然后，对每个采样比设定字典${\boldsymbol{\varPhi}}$，观测向量${\boldsymbol{y}}$由${\boldsymbol{\varPhi}} {\boldsymbol{\theta}}$与复高斯白噪声叠加获得，信噪比由0～15 dB，对平面的每个点，均进行25次蒙特卡罗实验，最终利用NMSE刻画信号重构准确率。

图 2  稀疏向量θ

下载: 全尺寸图片 幻灯片

图3为OMP, GD及VBI算法的相变图。其中各子图右上方区域代表可重建区域，左下方代表不可重建区域。可以看出，VBI算法相变图的可重建区域明显大于OMP和GD算法，因此该算法对信号的重构能力优于OMP和GD算法。

为进一步比较3种算法性能，以采样比为0.5观测稀疏向量${\boldsymbol{\theta}}$获取观测向量，并添加复高斯白噪声，信噪比由0～15 dB，并对每个信噪比设置了100个不同的噪声状态^[12]。用均方误差(Mean Square Error, MSE)衡量3种方法重构性能。3种方法不同信噪比下重构结果均方误差如图4所示。可以看出，3种算法重构误差均随信噪比增加而减小，与其他两种方法相比，VBI算法在各信噪比下的均方误差均为最小，表明VBI算法具有更好的重构性能与鲁棒性。此外OMP算法重构结果均方误差最大，GD算法次之。

图 3  3种算法相应相变图

下载: 全尺寸图片 幻灯片

图 4  3种算法平均重构误差随SNR变化曲线

下载: 全尺寸图片 幻灯片

表2提供了OMP, GD及VBI算法的运算复杂度对比，其中${k_0}$为稀疏度，$L$为信号长度。由表2可见，OMP算法运算复杂度最低，但低信噪比下重构误差最大，虽然GD算法与VBI算法运算复杂度相同，但VBI算法在运算复杂度可接受的情况下具有更优的重构效果。同时，可通过并行计算提高运算速度。

表 2  运算复杂度对比

重构算法	OMP	GD	VBI
运算复杂度	${\cal{O}}\left( {{k_0}{L^2}} \right)$	${\cal{O}}\left( {{L^{\rm{3}}}} \right)$	${\cal{O}}\left( {{L^{\rm{3}}}} \right)$

下载: 导出CSV 

| 显示表格

4.2   卫星数据仿真实验

本节通过仿真数据验证所提径向运动参数估计算法及高分辨成像算法的有效性。根据表3列出的参数生成卫星散射点模型(如图5所示)回波，假定目标剩余速度$\Delta {v_{\rm{R}}}$和剩余加速度$\Delta {a_{\rm{R}}}$分别为9 m/s和1 m/s²。稀疏步进调频信号共包含128个脉冲串，每个脉冲串包含从80个连续全频带脉冲中随机选择的64个脉冲(波形1)。通过对回波添加复高斯白噪声，回波信噪比以5 dB为步长从0 dB增加到15 dB，针对每个信噪比进行了25次不同噪声状态的独立试验，遗传算法种群数设置为40，遗传终止次数设置为20。将所提算法并与基于参数化字典的PSO算法^[10](算法1)进行对比。其中，算法1将距离像熵与平均距离像熵加权后作为目标函数。

表 3  雷达系统参数

${f_{\rm{c}}}$	${\rm{PRF}}$	${T_{\rm{R}}}$	$B$	$\Delta f$
10 GHz	6.4 kHz	20 μs	800 MHz	10 MHz

下载: 导出CSV 

| 显示表格

图 5  目标散射点分布图

下载: 全尺寸图片 幻灯片

目标剩余速度、加速度估计误差随SNR的变化如图6所示。可以看出，低信噪比下算法1的估计误差较大，不满足式(28a)及式(29)。这是由低信噪比条件下OMP算法较大的重构误差大导致的。随着信噪比提高，算法1估计误差有所下降，但估计结果仍不满足精度要求，会导致图像散焦。与算法1相比，所提算法在各信噪比下均具有鲁棒性，剩余速度与剩余加速度估计误差为1×10^–2 m/s和5×10^–3 m/s²，满足式(28a)及式(29)估计精度要求。低信噪比条件下，所提算法对运动参数估计的优越性能得益于：(1)基于Gamma-Gauss先验的重构算法能够实现HRRP的高精度重构；(2)遗传算法的全局优化能力。

图 6  速度加速度估计误差随SNR变化曲线

下载: 全尺寸图片 幻灯片

由于信噪比0 dB时算法1估计误差过大，无法实现聚焦成像，因此为公平比较，这里根据所提算法得到的运动参数估计值构造字典，进而分别采用OMP, GD, VBI算法求解式(10)，所得图像如图7所示。由图可见，所提算法对应成像结果轮廓清晰，虚假点少，聚焦效果均优于其他算法。特别地，所提算法可以更好地描述太阳能帆板细节。此外，3种方法对应的图像熵分别为0.3310, 0.2935, 0.2915。

图 7  仿真数据波形1成像结果

下载: 全尺寸图片 幻灯片

为了进一步验证所提算法的有效性，从80个连续全频带脉冲中随机选择48个脉冲生成新的稀疏步进调频信号(波形2)。同样，回波信噪比为0 dB。此时，所提算法的估计值$\left( {\Delta {{\hat v}_{\rm{R}}},\Delta {{\hat a}_{\rm{R}}}} \right)$为(8.9719, 0.9925)，算法1估计值为(8.1537, 0.4765)。由于算法1估计误差太大而无法聚焦成像，因此仅给出利用所提算法构造字典后，采用OMP, GD, VBI得到的成像结果。如图8所示，当脉冲数减少时，图像变得稀疏，但所提算法仍具有最好的成像效果，图像熵分别为0.3132, 0.2808, 0.2641。

图 8  仿真数据波形2成像结果

下载: 全尺寸图片 幻灯片

4.3   实测数据实验

本节通过Yak-42飞机实测数据验证所提算法的有效性。首先，Yak-42飞机实测数据为线性调频信号回波，通过对回波沿距离维等间隔分组，可以得到其步进调频信号回波；接着，在回波中加入由剩余运动引入的相位与复高斯白噪声；最终，通过对脉冲进行抽取获得稀疏步进调频信号。该数据的产生过程与本文算法的推导过程一致。图9为Yak-42飞机全频带实测数据成像结果，其中虚线圈代表弱散射点位置。

图 9  飞机全频带实测数据成像结果

下载: 全尺寸图片 幻灯片

对于平动补偿后的Yak-42飞机实测数据，分别将目标剩余速度和加速度设置为10 m/s和1 m/s²。回波全带宽为400 MHz，由64个连续脉冲信号形成步进调频信号，距离向采样点为256，方位向采样点为512。实验中，从全频带回波抽取出两种稀疏步进调频信号并记为波形1、波形2。波形1包含50个脉冲，波形2包含40个脉冲。通过加入复高斯白噪声使回波信噪比变为0 dB。首先，分别采用所提算法和算法1对回波进行运动参数估计，估计值如表4和表5所示，可见本文算法的估计值更接近目标真实参数。

表 4  Yak-42飞机剩余速度估计值表

	本文算法	算法1
波形1(m/s)	9.9763	9.8973
波形2(m/s)	9.9619	9.7503

下载: 导出CSV 

| 显示表格

表 5  Yak-42飞机剩余加速度估计值表

	本文算法	算法1
波形1(m/s2)	0.9975	1.1147
波形2(m/s2)	0.9994	1.2095

下载: 导出CSV 

| 显示表格

表 6  不同算法对Yak-42飞机实测数据成像的图像熵

	OMP	GD	VBI
波形1	0.2757	0.2108	0.1963
波形2	0.2086	0.1957	0.1767

下载: 导出CSV 

| 显示表格

然后，分别采用OMP, GD, VBI算法对波形1与波形2回波进行成像，结果如图10与图11所示。由图可见，所提算法可获得聚焦良好的图像，虚假点更少，并能够更好地保留飞机的轮廓和结构信息。与OMP算法相比，所提算法成像结果虚假点更少，与GD算法相比，所提算法可以保留更多的弱散射点。表6表明，两种波形下所提算法对应的图像熵均为最小。同时，由于波形2比波形1子脉冲更少，因此波形2的图像熵小于波形1的图像熵。

图 10  飞机实测数据波形1成像结果

下载: 全尺寸图片 幻灯片

图 11  飞机实测数据波形2成像结果

下载: 全尺寸图片 幻灯片

5.   结束语

本文针对稀疏步进调频信号对目标径向运动敏感且低信噪比下难以聚焦成像的问题，提出一种低信噪比稀疏步进调频信号运动参数估计与高分辨成像算法。针对HRRP合成，引入Gamma-Gauss共轭先验构建概率图模型，进而采用VBI算法对其进行推断求解。由于充分利用了散射点和噪声的统计特性，该方法噪声鲁棒性，HRRP重构精度高。同时，构造参数化字典，并以图像熵最小为准则，采用遗传算法实现对运动参数和HRRP的准确估计。未来工作将重点研究基于环境感知的稀疏步进调频信号波形设计与复杂运动目标高分辨成像方法。