CN113534151B - 基于离网稀疏贝叶斯学习的双频段isar成像方法 - Google Patents

Abstract

本发明属于雷达信号处理技术领域，公开了一种基于离网稀疏贝叶斯学习的双频段ISAR成像方法，包括：对高低频段回波数据进行预处理，得到高低频段观测数据矩阵；根据高低频段观测数据矩阵建立高低频段数据的离网稀疏表征模型，并基于二阶泰勒展开简化高低频段数据的离网稀疏表征模型；根据高低频段在全频段的索引向量，生成对应的基字典矩阵和离网补偿字典矩阵；利用稀疏贝叶斯学习算法依次对不同方位单元中所包含的离网散射中心位置和散射强度进行求解，即可实现目标的二维高分辨成像。解决了观测目标上等效散射中心偏离字典网格时现有成像算法旁瓣高难以聚焦的问题，且在低信噪比下仍可得到聚焦的合成高分辨二维像，具有一定的噪声稳健性。

Description

基于离网稀疏贝叶斯学习的双频段ISAR成像方法

技术领域

本发明涉及雷达信号处理技术领域，具体涉及基于离网稀疏贝叶斯学习的双频段ISAR成像方法。

背景技术

逆合成孔径雷达(Inverse Synthetic Aperture Radar，ISAR)凭借其全天时、全天候、远距离和高分辨成像的能力，广泛应用于目标探测和态势感知等军事与民用领域。ISAR图像距离维分辨率与雷达发射信号的带宽有关，带宽越大，分辨率越高。因此，为了获得高的距离分辨率，通常需要雷达具备大带宽信号的发射和接收处理能力。然而雷达发射信号的带宽并不是无限增加的，而是受到雷达收发链路硬件设备性能的制约。而多频带相干合成技术是解决上述问题的有效途径，在不改变雷达现有硬件的情况下，基于不同频段雷达对同一个目标相同观测视角的回波，通过多频段相干合成技术，即可实现ISAR成像距离分辨率的显著提升。

目前多频带相干合成ISAR成像方法主要分为两大类。一类是基于传统谱估计的多频段相干合成ISAR成像方法，该类方法在低信噪比下难以获得较高的精度；第二类是基于稀疏表示的多频段相干合成ISAR成像方法。该类方法又可细分为三种：第一种是基于贪婪算法的多频段相干合成成像方法，如最常见的正交匹配追踪(Orthogonal MatchingPursuit，OMP)算法，这类算法优点是速度快，但精度不高，且易受噪声影响。第二种是基于L_p范数优化的多频段相干合成ISAR成像方法，如较典型的L₁范数正则化方法、平滑L₀算法等；该类方法虽然求解精度高于贪婪类算法，但是需要手动调节正则化项参数，算法性能依赖于手动调参。第三种是基于稀疏贝叶斯学习(Sparse Bayesian Learning，SBL)的多频段相干合成ISAR成像方法，也是目前多频带相干合成ISAR成像方法的研究热点；该类方法不需要手动调节参数，而是精确建模散射中心和噪声的先验分布，通过相应的迭代求解框架自适应求解散射中心位置与散射强度，从而得到更精确的估计结果。

现有稀疏贝叶斯学习方法一般通过构建单层字典进行双频段合成ISAR成像，而在实际ISAR成像场景中，目标上的等效散射中心很可能偏离该单层字典划分的网格，从而导致成像结果出现高旁瓣和虚假的散射中心，使得成像结果散焦和分辨率下降。针对上述问题，有学者通过构建多级动态字典划分更细的网格，但是该方法求得的稀疏向量各分量并非处在均匀的网格上，所以无法直接利用该稀疏向量作为目标一维距离像，而是需要恢复全频段回波数据，然后再进行距离向傅里叶变换得到均匀网格下的一维距离像，仍然无法避免散射中心位置无法落在傅里叶网格上的问题，导致成像结果散焦。

发明内容

为了克服现有多频段相干合成成像方法在低信噪比和等效散射中心偏离网格时无法准确获得目标高分辨二维聚焦图像的问题，本发明提出了基于离网稀疏贝叶斯学习的双频段ISAR成像方法，解决了观测目标上等效散射中心偏离字典网格时现有成像算法旁瓣高难以聚焦的问题，并且，本发明在低信噪比下仍可得到聚焦的合成高分辨二维像，具有一定的噪声稳健性。

为了达到上述目的，本发明采用以下技术方案予以实现。

基于离网稀疏贝叶斯学习的双频段ISAR成像方法，包括以下步骤：

步骤1，对高低频段回波数据进行预处理，得到高低频段观测数据矩阵；

步骤2，根据所述高低频段观测数据矩阵建立高低频段数据的离网稀疏表征模型，并基于二阶泰勒展开简化所述高低频段数据的离网稀疏表征模型，得到基于二阶泰勒展开的双频段离网稀疏表征模型；根据高低频段在全频段的索引向量，生成对应的基字典矩阵和离网补偿字典矩阵；

步骤3，利用稀疏贝叶斯学习算法依次对不同方位单元中所包含的离网散射中心位置和散射强度进行求解，即可实现目标的二维高分辨成像。

本发明技术方案的特点和进一步的改进为：

(1)步骤1具体包含以下子步骤：

子步骤1.1，采用解线频调方式对高低频段回波数据进行脉冲压缩处理，再对脉冲压缩后的高低频段回波数据进行平动补偿，得到高低频段距离时域的回波

为：

其中，

是距离快时间，t_m是方位慢时间，K为观测目标上等效散射中心的个数，σ_k为第k个散射中心的后向散射系数，R_k(t_m)为第m次回波第k个散射中心与雷达之间的距离，R_ref为参考距离，ΔR_k(t_m)＝R_k(t_m)-R_ref，C为电磁波在真空中的传播速度，T_p为发射信号脉冲宽度，γ为发射信号调频率，f_ci为发射信号载频，j是复数单位；

子步骤1.2，经过脉冲压缩处理后，不同散射中心的时域回波变为单频脉冲，其频率为

并且/>

其中x_k、y_k和w分别表示散射中心的横坐标、纵坐标与等效转动角速度；将上式进行转换，将所述高低频段距离时域回波表达式写为：

其中，f_r1和f_r2分别表示低频段和高频段雷达距离时域回波信号对应的距离频率，f_c1和f_c2分别表示低频段和高频段雷达距离时域回波信号的载频，B₁和B₂分别表示低频段和高频段雷达距离时域回波信号对应的带宽；

子步骤1.3，将高低频段回波合成之后的回波称为合成全频段回波，利用合成全频段回波对应的中心频率f_c，对所述高低频段距离时域回波s(f_ri，t_m)进行统一的Keystone变换，消除越距离单元徙动和高低频段观测数据载频不同引起的方位分辨率差异，再进行方位向脉压，得到目标高低频段方位成像结果s(f_ri，f_d)：

其中，f_d为方位频率；

子步骤1.4，将所述目标高低频段方位成像结果s(f_ri，f_d)沿方位向进行截取，得到包含目标的距离单元；再沿距离数据域对所述包含目标的距离单元进行截取，得到距离向截取后的高低频段回波数据s(f_ri，f_d)：

其中，

f_d取值为方位向目标所在区域对应的频率，距离向的截取操作使得/>

利用公式/>

得到子频段截取后的距离向点数，其中N_Ri为子频段截取前的距离向采样点数，f_si为第i个频段的采样频率；

子步骤1.5，高低频段的频率采样间隔Δf_i与其对应的带宽B_i、脉冲宽度T_pi、采样频率f_si有关，即

对距离向截取后的高低频段回波数据进行插值、抽取处理调节频率采样间隔，得到距离向网格分辨率一致的高低频段数据；其中，高低频段距离向点数变为N_i，i＝1，2，Δf₁＝Δf₂，距离网格分辨率统一为Δf，记f_ri+f_ci＝f_c+n_iΔf，其中n_i为子频带i在合成全频段中对应的索引；记方位采样索引为m，采样间隔为Δf_d，则预处理后的高低频段回波数据为：

其中，

该系数与i无关，于是高低频段的回波数据已被校正到统一的回波模型下；

子步骤1.6，根据预处理后的高低频段回波数据构建高低频段观测数据矩阵Y＝[Y_L；Y_H]，其维度为N×M，其中Y_L＝[s(n₁，m₁)，…，s(n₁，m_M)]为低频段回波数据，Y_H＝[s(n₂，m₁)，…，s(n₂，m_M)]为高频段回波数据，N＝N₁+N₂为高低频段数据距离点数之和，M为方位截取后的回波次数。

(2)步骤2包含以下子步骤：

子步骤2.1，若使用相位项2πq表示任意一个处于[-π，π)的相位，根据指数函数的周期性可知，对于目标的每一个等效散射中心k，存在q_k使得

即q包含了目标等效散射中心位置的所有可能取值；其中/>

根据稀疏表示理论，将q等分为D份，用

表示第k个散射中心对应的真实位置，用d_k表示离真实位置/>

最近的网格点，用δ_k表示真实位置与其最近的网格点之间的距离，则

从而构建出高低频段数据的离网稀疏表征模型：

Y＝ΦX+E

其中，E为噪声矩阵；稀疏矩阵X为高低频段回波相干合成的二维成像结果，X的维度为D×M；离网字典

其维度为N×D；/>

为目标真实位置/>

的所有可能取值；离网字典Φ的每一个原子表示为：

其中，n为高低频段在合成全频段中的索引向量，其维度为N×1；n₁为低频段在合成全频段中的索引向量；n₂为高频段在合成全频段中的索引向量；

其中

稀疏矩阵X中的每一列x都是稀疏向量，其不为零项的值即表示目标散射中心的幅度σ″，不为零项对应的原子相位即对应目标散射中心的位置y_k，于是求解出的稀疏向量x即可表示目标一维距离像，而所有回波求解出的稀疏矩阵X即为高低频段回波相干合成的二维成像结果；

子步骤2.2，由于δ_r取值很小，将原子

在d_r处进行二阶泰勒展开，即

其中

若记基字典矩阵为A＝[a(d₁)，…，a(d_D)]，一阶补偿字典矩阵B＝[b(d₁)，…，b(d_D)]，二阶补偿字典矩阵C＝[c(d₁)，…，c(d_D)]，其中δ＝[δ₁，…，δ_D]^T称为模型补偿向量，于是得到基于二阶泰勒展开的双频段离网稀疏表征模型：

Y＝ΦX+E＝(A+Bdiag(δ)+Cdiag(δ²))X+E

其中

表示两个向量的hadamard积，diag(δ)表示以δ为对角线元素的对角矩阵；根据高低频段在全频段的索引向量，生成对应的基字典矩阵A、一阶补偿字典矩阵B、二阶补偿字典矩阵C。

(3)步骤3包含以下子步骤：

子步骤3.1，第m次回波的基于二阶泰勒展开的双频段离网稀疏表征模型为：

y＝Φx+ε＝(A+Bdiag(δ)+Cdiag(δ²))x+ε

其中，观测向量y为高低频段观测数据矩阵Y的某一列，稀疏向量x为稀疏矩阵X的某一列；设噪声向量ε的实部和虚部均独立同分布于均值为0、协方差为

的高斯分布，其中β为噪声精度，I为单位矩阵，于是ε就服从于均值为0、协方差为β^-1I的复高斯分布，可以得到已知β的条件下ε的概率密度函数p(ε|β)为：

p(ε|β)＝CN(ε|0，β^-1I)

其中，CN()表示复高斯分布的概率密度函数，设噪声精度β服从于参数为(c，d)的伽马分布，于是得到β的概率密度函数p(β)为：

p(β)＝Gam(β|c，d)

其中，Gam()表示伽马分布的概率密度函数，c、d分别表示伽马分布的形状参数和逆尺度参数；

观测向量y的似然概率p(y|x，β，δ)服从均值为Φx、方差为β^-1I的复高斯分布，即：

p(y|x，β，δ)＝CN(Φx，β^-1I)

将散射中心位置和幅度向量建模为两层的高斯逆伽马先验，即认为x服从于方差为α的复高斯分布，于是得到已知α的条件下x的概率密度函数p(x|α)为：

p(x|α)＝CN(0，Λ^-1I)

其中，方差矩阵Λ＝diag(α)，α＝[α₁，…，α_r，…，α_D]^T∈R^D，并且α_r独立同分布于参数为(a，b)的逆伽马分布，即α_r的概率密度函数p(α_r)为：

p(α_r)＝IGam(α_r|a，b)，(r＝1，...，D)

其中，IGam()表示逆伽马分布的概率密度函数，a、b分别表示逆伽马分布的形状参数和尺度参数；假设补偿向量δ_r服从于均匀分布，即δ_r的概率密度函数p(δ_r)表示为：

子步骤3.2，利用期望最大化EM框架对未知参数进行贝叶斯推断，基于证据最大化理论估计参数θ＝(α，β，δ)，估计值

表示为：

其中，p(θ|y)为已知y的条件下θ的概率密度函数，p(y，θ)为y，θ的联合概率密度函数，即最大化(y，α，β，δ)的联合概率密度函数的对数，其表达式为：

其中，D＝β^-1I+ΦΛΦ为边缘似然分布的协方差矩阵；

子步骤3.3，利用期望最大化EM框架通过最大化代价函数的下界Q(α，β，δ)求参数，即

maxQ(α，β，δ)＝max E_{x～p(x|y，α，β，δ)}{ln p(y，x，α，β，δ)}

下界Q(α，β，δ)等于完全数据的对数联合分布函数ln p(y，x，α，β，δ)对隐变量x的条件概率分布p(x|y，α，β，δ)的期望；根据概率图模型可写出其联合分布为：

p(y，x，α，β，δ)＝p(y|x，β，δ)p(x|α)p(α)p(β)p(δ)

其中，p(δ)为δ的概率密度函数；

E步：求隐变量x的后验期望，x的后验概率分布为

p(x|y，α，β，δ)∝p(y|x，β，δ)p(x|α)＝CN(x|μ，∑)

将期望值μ作为x的估计值，∑是x的后验协方差矩阵；

M步：通过最大化E{log p(x|α)p(α)}得到方差α的更新公式，通过最大化E{log p(y|x，β，δ)p(β)}得到噪声精度β的更新公式，通过最大化E{log p(y|x，β，δ)p(δ)}得到补偿因子δ的更新公式；

子步骤3.4，利用所述期望最大化EM框架依次对不同方位单元中所包含的散射中心位置和散射强度进行求解，即可得到目标的二维ISAR图像。

(4)子步骤3.4包括以下子步骤：

子步骤3.4.1，初始化超参数a＝b＝10^-4，模型补偿向量δ＝{0}_D×1，模型稀疏度为K_m，最大迭代次数Iter_max＝3000，迭代收敛门限τ＝10^-2；

子步骤3.4.2，对于高低频段观测数据矩阵的第m列，初始化噪声精度

其中/>

S_i为子频段i未截取前的二维像，n_o∈I_noise表示噪声所在的距离单元索引，N_i为高低频段距离向单元数，并初始化为c＝1，d＝Nβ_m；初始化稀疏向量x的方差α＝|A^Hy|/N；当前迭代次数iter＝1；

子步骤3.4.3，根据如下公式更新稀疏向量x：

μ＝β∑Φ^Hy

∑＝(βΦ^HΦ+Λ^-1)^-1

将期望值μ作为x的估计值；

子步骤3.4.4，根据如下公式更新α_r：

其中，<|x_r|²〉＝E{|x_r|²}＝|μ_r|²+Re(∑_rr)，μ_r为μ的第r个元素，∑_rr为∑的第r行第r列元素，||表示取模值操作，Re()表示取实部操作；

根据如下公式更新噪声精度β：

其中，

子步骤3.4.5，找到前K_m个最大的α对应的索引I_k，并且通过求方程f′(δ_k)(k∈I_k)的根更新每一个δ_k；其中，f′(δ_k)通过最大化E{log p(y|x，β，δ)p(δ)}推导出来的：

其中

v₁＝Re[diag(μ^*)B^H(y-Aμ)]-Re[diag(B^HA∑)]

v₂＝Re[diag(μ^*)C^H(y-Aμ)]-Re[diag(C^HA∑)]

()^*表示取共轭操作，()^H表示取共轭转置操作，()_k，k表示取出矩阵的第k行第k列元素的操作，

表示取出矩阵的第k行除掉第k列外其他元素的操作，/>

表示取出矩阵的第k列除掉第k行其他元素的操作，()^-k表示取出向量除第k个元素外其余元素的操作，()_k表示取出向量第k个元素的操作；

子步骤3.4.6，判断迭代终止条件

或者iter＜Iter_max是否满足，若有一个条件满足则停止迭代；否则，令迭代次数iter增加1，转到子步骤3.4.3；

子步骤3.4.7，迭代完成后，利用模型量化误差求解目标离网位置

其中/>

表示索引I_k对应的补偿量；然后构造离网字典

利用子步骤3.4.3求得的μ作为第m次回波的高分辨一维距离像；

子步骤3.4.8，判断当前方位单元序号是否等于M，若是，输出高分辨二维像；否则，执行子步骤3.4.2。

与现有技术相比，本发明的有益效果为：

本发明的离网稀疏贝叶斯学习算法首先基于二阶泰勒展开构建离网稀疏表示模型，然后引入高斯-逆伽马先验作为稀疏向量和噪声的先验，接着采用EM框架进行后验推断，最后将求得的稀疏向量作为目标的一维距离像。该发明克服了现有算法针对离网目标成像出现散焦的问题，并且无需恢复全频段回波数据，克服了传统通过傅里叶变换成像的缺点，在低信噪比环境下能够获得精确的高分辨二维像。

附图说明

下面结合附图和具体实施例对本发明做进一步详细说明。

图1为本发明实施例提供的基于离网稀疏贝叶斯学习的双频段ISAR成像方法流程图；

图2为基于二阶泰勒展开模型的离网稀疏贝叶斯学习算法流程图；

图3为飞机二维点目标模型图；

图4为本发明的仿真结果图；其中图4(a)为低频段二维成像结果图；图4(b)为高频段二维成像结果图；图4(c)为真实全频段二维成像结果图；

图4(d)为采用本发明重构的二维成像结果图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1为本发明实施例提供的基于离网稀疏贝叶斯学习算法的双频段合成ISAR成像方法的流程示意图，参考图1，本发明的基于离网稀疏贝叶斯学习的双频段ISAR成像方法，包括以下步骤：

步骤1，对高低频段回波数据进行预处理，得到高低频段观测数据矩阵。

具体的，步骤1包含以下子步骤：

子步骤1.1，采用解线频调方式对高低频段回波数据进行脉冲压缩处理，将不同散射中心对应的单频脉冲在时域上对齐，并且对脉冲压缩处理后的高低频段一维距离像数据进行平动补偿，校正高低频段回波中的包络偏移和初相误差，可得到目标距离时域的回波如下：

其中，

是距离快时间，t_m是方位慢时间，K为观测目标上等效散射中心的个数，σ_k为第k个散射中心的后向散射系数，R_k(t_m)为第m次回波第k个散射中心与雷达之间的距离，R_ref是参考距离，ΔR_k(t_m)＝R_k(t_m)-R_ref，C为电磁波在真空中的传播速度，T_p为发射信号脉冲宽度，γ为发射信号调频率，f_ci为发射信号载频，j是复数单位。

子步骤1.2，经过解线频调脉冲压缩处理后，不同散射中心的时域回波变为单频脉冲，其频率为

并且/>

其中x_k、y_k和w分别表示散射中心的横坐标、纵坐标与等效转动角速度。将上式进行转换，可以将高低频段距离时域回波表达式写为：

其中，f_r1和f_r2分别表示低频段和高频段雷达回波信号对应的距离频率，f_c1和f_c2分别表示低频段和高频段雷达回波信号的载频，B₁和B₂分别表示低频段和高频段雷达回波信号对应的带宽，上式所示回波称为距离数据域回波。

子步骤1.3，将高低频段回波合成之后的回波称为合成全频段回波，利用合成全频段回波对应的中心频率f_c，对高低频段距离时域回波进行统一的Keystone变换，消除越距离单元徙动和高低频段观测数据载频不同引起的方位分辨率差异，得到方位向分辨率一致的高低频段数据。之后再进行方位向脉压得到目标高低频段方位成像结果：

其中，f_d为方位频率，f_c为合成全频段回波的中心频率，sinc函数表达式为sinc(x)＝sin(x)/x，rect函数表达式为

子步骤1.4，将目标高低频段方位成像结果沿方位向进行截取，保留观测目标所在区域，降低后续相干合成成像计算量。然后沿距离数据域对包含目标的距离单元进行截取操作，以便得到有效的数据进行后续的合成操作，此时高低频段数据可表示为：

其中

f_d取值为方位向目标所在区域对应的频率，距离向的截取操作使得/>

利用公式/>

得到子频段截取后的距离向点数，其中N_Ri为子频段截取前的距离向采样点数，f_si为第i个频段的采样频率。

子步骤1.5，高低频段的频率采样间隔Δf_i与其对应的带宽B_i、脉冲宽度T_pi、采样频率f_si有关，即

若要进行高低频段的合成，需保证距离网格分辨率Δf₁＝Δf₂。所以当高低频段信号的带宽、脉冲宽度、采样频率不同导致距离网格分辨率不同时，需要对高低频段数据进行插值、抽取处理调节采样间隔，以得到距离向网格分辨率一致的高低频段数据。此时，高低频段距离向点数变为N_i，i＝1，2，距离网格分辨率统一为Δf。记f_ri+f_ci＝f_c+n_iΔf，其中n_i为子频带i在合成全频段中对应的索引。记方位采样索引为m，采样间隔为Δf_d，于是单次回波模型可简化为：

其中

该系数与i无关，于是高低频段的回波数据已被校正到统一的回波模型下。

子步骤1.6，根据预处理后的高低频段回波数据构建高低频段观测数据矩阵Y＝[Y_L；Y_H]，其维度为N×M，其中Y_L＝[s(n₁，m₁)，…，s(n₁，m_M)]为低频段回波数据，Y_H＝[s(n₂，m₁)，…，s(n₂，m_M)]为高频段回波数据，N＝N₁+N₂为高低频段数据距离点数之和，M为方位截取后的回波次数。

步骤2，根据所述高低频段观测数据矩阵建立高低频段数据的离网稀疏表征模型，并基于二阶泰勒展开简化所述高低频段数据的离网稀疏表征模型，得到基于二阶泰勒展开的双频段离网稀疏表征模型；根据高低频段在全频段的索引向量，生成对应的基字典矩阵和离网补偿字典矩阵。

具体的，步骤2包含以下子步骤：

子步骤2.1，若使用相位项2πq(其中

)表示任意一个处于[-π，π)的相位，根据指数函数的周期性可知，对于目标的每一个等效散射中心k，存在q_k使得

即q包含了目标等效散射中心位置的所有可能取值。根据稀疏表示理论，本发明将q等分为D份，在网稀疏表示模型认为目标等效散射中心均处于预先划分的网格上，但是实际中目标等效散射中心可能偏离网格。如果用/>

表示第k个散射中心对应的真实位置，用d_k表示离真实位置/>

最近的网格点，用δ_k表示真实位置与其最近的网格点之间的距离，于是/>

从而构建出高低频段数据的离网稀疏表征模型：

Y＝ΦX+E

其中，E为噪声矩阵；稀疏矩阵X为高低频段回波相干合成的二维成像结果，X的维度为D×M。

称为离网字典，其维度为N×D，/>

为目标真实位置/>

的所有可能取值，离网字典Φ的每一个原子表示为：

其中，n为高低频段在合成全频段中的索引向量，其维度为N×1；n₁为低频段在合成全频段中的索引向量；n₂为高频段在合成全频段中的索引向量；

其中

X的维度为D×M，每一列x都是稀疏向量，其不为零项的值即表示目标散射中心的幅度σ″，不为零项对应的原子相位即对应目标散射中心的位置y_k，于是求解出的稀疏向量x即可表示目标一维距离像，而所有回波求解出的稀疏矩阵X即为高低频段回波相干合成的二维成像结果。

子步骤2.2，由于δ_r取值很小，可以将原子

在d_r处进行二阶泰勒展开，即

其中

若记基字典矩阵为A＝[a(d₁)，…，a(d_D)]，一阶补偿字典矩阵B＝[b(d₁)，…，b(d_D)]，二阶补偿字典矩阵C＝[c(d₁)，…，c(d_D)]，其中δ＝[δ₁，…，δ_D]^T称为模型补偿向量，于是得到基于二阶泰勒展开的双频段离网稀疏表征模型：

Y＝ΦX+E＝(A+Bdiag(δ)+Cdiag(δ²))X+E

其中

表示两个向量的hadamard积，diag(δ)表示以δ为对角线元素的对角矩阵。根据高低频段在全频段的索引向量，生成对应的基字典矩阵A、一阶补偿字典矩阵B、二阶补偿字典矩阵C。

步骤3，利用稀疏贝叶斯学习算法依次对不同方位单元中所包含的离网散射中心位置和散射强度进行求解，即可实现目标的二维高分辨成像。

步骤3包含以下子步骤：

子步骤3.1，首先针对某次回波的基于二阶泰勒展开的双频段离网稀疏表征模型建立对应的概率图模型。该次回波的基于二阶泰勒展开的双频段离网稀疏表征模型为

y＝Φx+ε＝(A+Bdiag(δ)+Cdiag(δ²))x+ε

其中y为Y的某一列，x为X的某一列。假设噪声向量ε的实部和虚部均独立同分布于均值为0、协方差为

的高斯分布，其中β为噪声精度，I为单位矩阵，于是ε就服从于均值为0、协方差为β^-1I的复高斯分布，可以得到已知β的条件下ε的概率密度函数p(ε|β)为：

p(ε|β)＝CN(ε|0，β^-1I)

CN()表示复高斯分布的概率密度函数，并且假设β服从于参数为(c，d)的伽马分布，于是得到β的概率密度函数p(β)为：

p(β)＝Gam(β|c，d)

其中Gam()表示伽马分布的概率密度函数，c，d分别表示伽马分布的形状参数和逆尺度参数。并且，该噪声先验使得本发明所提离网稀疏贝叶斯学习算法对噪声具有一定的鲁棒性。

观测向量y的似然概率p(y|x，β，δ)服从均值为Φx、方差为β^-1I的复高斯分布，即：

p(y|x，β，δ)＝CN(Φx，β^-1I)

将散射中心位置和幅度向量建模为两层的高斯逆伽马先验，即认为x服从于方差为α的复高斯分布，于是得到已知α的条件下x的概率密度函数p(x|α)为：

p(x|α)＝CN(0，Λ^-1I)

其中，方差矩阵Λ＝diag(α)，α＝[α₁，…，α_r，…，α_D]^T∈R^D，并且α_r(r＝1，...，D)独立同分布于参数为(a，b)的逆伽马分布，即α_r的概率密度函数p(α_r)为：

p(α_r)＝IGam(α_r|a，b)

IGam()表示逆伽马分布的概率密度函数，a，b分别表示逆伽马分布的形状参数和尺度参数。假设补偿向量δ_r(r＝1，...，D)服从于均匀分布，即δ_r的概率密度函数p(δ_r)表示为：

子步骤3.2，利用期望最大化(Expectation-Maximization，EM)框架对未知参数进行贝叶斯推断。本发明基于证据最大化理论估计参数θ＝(α，β，δ)，估计值

表示为：

其中，p(θ|y)为已知y的条件下θ的概率密度函数，p(y，θ)为y，θ的联合概率密度函数，即最大化(y，α，β，δ)的联合概率密度函数的对数，其表达式为

其中，D＝β^-1I+ΦΛΦ为边缘似然分布的协方差矩阵。

子步骤3.3，由于直接通过最大化代价函数lnp(y，α，β，δ)难以得到参数的闭式解，所以采用基于期望最大化(Expectation-Maximization，EM)的迭代算法估计参数。EM算法用于含有隐变量的概率模型参数的极大似然估计或极大后验概率估计，该算法通过最大化代价函数的下界Q(α，β，δ)求参数，即

max Q(α，β，δ)＝max E_{x～p(x|y，α，β，δ)}{ln p(y，x，α，β，δ)}

下界Q(α，β，δ)等于完全数据的对数联合分布函数ln p(y，x，α，β，δ)对隐变量x的条件概率分布p(x|y，α，β，δ)的期望。根据概率图模型可写出其联合分布为

p(y，x，α，β，δ)＝p(y|x，β，δ)p(x|α)p(α)p(β)p(δ)

其中，p(δ)为δ的概率密度函数；

E步：求隐变量x的后验期望，x的后验概率分布为

p(x|y，α，β，δ)∝p(y|x，β，δ)p(x|α)＝CN(x|μ，∑)

将期望值μ作为x的估计值，∑是x的后验协方差矩阵。

M步：通过最大化E{log p(x|α)p(α)}得到方差α的更新公式，通过最大化E{log p(y|x，β，δ)p(β)}得到噪声精度β的更新公式，通过最大化E{log p(y|x，β，δ)p(δ)}得到补偿因子δ的更新公式。

子步骤3.4，利用上述框架依次对不同方位单元中所包含的散射中心位置和散射强度进行求解，即可得到目标的二维ISAR图像。

图2为基于二阶泰勒展开模型的离网稀疏贝叶斯学习算法流程示意图；参照图2，子步骤3.4包括以下子步骤：

子步骤3.4.1，初始化超参数a＝b＝10^-4，模型补偿向量δ＝{0}_D×1，模型稀疏度为K_m，最大迭代次数Iter_max＝3000，迭代收敛门限τ＝10^-2。

子步骤3.4.2，对于高低频段观测数据矩阵的第m列，初始化噪声精度

其中/>

S_i为子频段i未截取前的二维像，n_o∈I_noise表示噪声所在的距离单元索引，Ni为高低频段距离向单元数，并初始化为c＝1，d＝Nβ_m；初始化稀疏向量x的方差为α＝|A^Hy|/N；当前迭代次数iter＝1。

子步骤3.4.3，根据如下公式更新稀疏向量x：

μ＝β∑Φ^Hy

∑＝(βΦ^HΦ+Λ^-1)^-1

将期望值μ作为x的估计值。

子步骤3.4.4，根据如下公式更新α_r：

其中〈|X_r|²>＝E{|x_r|²}＝|μ_r|²+Re(∑_rr)，μ_r为μ的第r个元素，∑_rr为∑的第r行第r列元素，||表示取模值操作，Re()表示取实部操作。

根据如下公式更新噪声精度β：

其中

子步骤3.4.5，找到前K_m个最大的α对应的索引I_k，并且通过求方程f′(δ_k)(k∈I_k)的根更新每一个δ_k。其中，f′(δ_k)通过最大化E{log p(y|x，β，δ)p(δ)}推导出来的。

其中

v₁＝Re[diag(μ^*)B^H(y-Aμ)]-Re[diag(B^HA∑)]

v₂＝Re[diag(μ^*)C^H(y-Aμ)]-Re[diag(C^HA∑)]

()^*表示取共轭操作，()^H表示取共轭转置操作。()_k，k表示取出矩阵的第k行第k列元素的操作，

表示取出矩阵的第k行除掉第k列外其他元素的操作，/>

表示取出矩阵的第k列除掉第k行其他元素的操作，()^-k表示取出向量除第k个元素外其余元素的操作，()_k表示取出向量第k个元素的操作。

子步骤3.4.6，判断迭代终止条件

或者iter＜Iter_max是否满足，若有一个条件满足则停止迭代；否则，令迭代次数iter增加1，即iter＝iter+1，转到子步骤3.4.3。/>

子步骤3.4.7，迭代完成后，利用模型量化误差求解目标离网位置

其中/>

表示索引I_k对应的补偿量。然后构造离网字典

利用子步骤3.4.3求得的μ作为第m次回波的高分辨一维距离像。

子步骤3.4.8，判断当前方位单元序号是否等于M，若是，输出高分辨二维像；否则，执行子步骤3.4.2。

进一步的，以下通过仿真实验对本发明实施例方法的效果进行验证：

1、仿真条件

飞机二维点目标模型如图3所示。仿真假设雷达发射线性调频信号，低频段雷达载频为10.5GHz，带宽为1GHz，高频段雷达载频为14GHz，带宽为2GHz，脉冲宽度均为500μs，采样频率均为1MHz，脉冲重复频率为500Hz。针对高低频段方位像数据加入信噪比为5dB的高斯噪声。

2、仿真实验内容及结果分析

根据图4(a)和图4(b)可以知道，单个频段由于带宽受限难以区分距离较近的目标散射中心，导致成像结果分辨率较低。图4(c)为真实全频段二维成像结果图，图4(d)为采用本发明重构的二维成像结果图，可以看到本发明所提算法能够有效抑制伪峰，并且能够克服傅里叶变换离散网格的限制，在低信噪比下仍能获得高分辨的二维ISAR图像。

本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

虽然，本说明书中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。

Claims (2)

1.基于离网稀疏贝叶斯学习的双频段ISAR成像方法，其特征在于，包括以下步骤：

步骤1，对高低频段回波数据进行预处理，得到高低频段观测数据矩阵；

步骤1具体包含以下子步骤：

子步骤1.1，采用解线频调方式对高低频段回波数据进行脉冲压缩处理，再对脉冲压缩后的高低频段回波数据进行平动补偿，得到高低频段距离时域的回波

为：

其中，

是距离快时间，t_m是方位慢时间，K为观测目标上等效散射中心的个数，σ_k为第k个散射中心的后向散射系数，R_k(t_m)为第m次回波第k个散射中心与雷达之间的距离，R_ref为参考距离，ΔR_k(t_m)＝R_k(t_m)-R_ref，C为电磁波在真空中的传播速度，T_p为发射信号脉冲宽度，γ为发射信号调频率，f_ci为发射信号载频，j是复数单位；

子步骤1.2，经过脉冲压缩处理后，不同散射中心的时域回波变为单频脉冲，其频率为

并且/>

其中x_k、y_k和w分别表示散射中心的横坐标、纵坐标与等效转动角速度；将上式进行转换，将所述高低频段距离时域回波表达式写为：

其中，f_r1和f_r2分别表示低频段和高频段雷达距离时域回波信号对应的距离频率，f_c1和f_c2分别表示低频段和高频段雷达距离时域回波信号的载频，B₁和B₂分别表示低频段和高频段雷达距离时域回波信号对应的带宽；

子步骤1.3，将高低频段回波合成之后的回波称为合成全频段回波，利用合成全频段回波对应的中心频率f_c，对所述高低频段距离时域回波s(f_ri，t_m)进行统一的Keystone变换，消除越距离单元徙动和高低频段观测数据载频不同引起的方位分辨率差异，再进行方位向脉压，得到目标高低频段方位成像结果s(f_ri，f_d)：

其中，f_d为方位频率；

子步骤1.4，将所述目标高低频段方位成像结果s(f_ri，f_d)沿方位向进行截取，得到包含目标的距离单元；再沿距离数据域对所述包含目标的距离单元进行截取，得到距离向截取后的高低频段回波数据s(f_ri，f_d)：

其中，

f_d取值为方位向目标所在区域对应的频率，距离向的截取操作使得/>

利用公式/>

得到子频段截取后的距离向点数，其中N_Ri为子频段截取前的距离向采样点数，f_si为第i个频段的采样频率；

子步骤1.5，高低频段的频率采样间隔Δf_i与其对应的带宽B_i、脉冲宽度T_pi、采样频率f_si有关，即

对距离向截取后的高低频段回波数据进行插值、抽取处理调节频率采样间隔，得到距离向网格分辨率一致的高低频段数据；其中，高低频段距离向点数变为N_i，i＝1，2，Δf₁＝Δf₂，距离网格分辨率统一为Δf，记f_ri+f_ci＝f_c+n_iΔf，其中n_i为子频带i在合成全频段中对应的索引；记方位采样索引为m，采样间隔为Δf_d，则预处理后的高低频段回波数据为：

其中，

该系数与i无关，于是高低频段的回波数据已被校正到统一的回波模型下；

子步骤1.6，根据预处理后的高低频段回波数据构建高低频段观测数据矩阵Y＝[Y_L；Y_H]，其维度为N×M，其中Y_L＝[s(n₁，m₁)，…，s(n₁，m_M)]为低频段回波数据，Y_H＝[s(n₂，m₁)，…，s(n₂，m_M)]为高频段回波数据，N＝N₁+N₂为高低频段数据距离点数之和，M为方位截取后的回波次数；

步骤2，根据所述高低频段观测数据矩阵建立高低频段数据的离网稀疏表征模型，并基于二阶泰勒展开简化所述高低频段数据的离网稀疏表征模型，得到基于二阶泰勒展开的双频段离网稀疏表征模型；根据高低频段在全频段的索引向量，生成对应的基字典矩阵和离网补偿字典矩阵；

步骤2包含以下子步骤：

子步骤2.1，若使用相位项2πq表示任意一个处于[-π，π)的相位，根据指数函数的周期性可知，对于目标的每一个等效散射中心k，存在q_k使得

即q包含了目标等效散射中心位置的所有可能取值；其中/>

C为电磁波在真空中的传播速度，y_k为散射中心的纵坐标，Δf为距离网格分辨率；

根据稀疏表示理论，将q等分为D份，用

表示第k个散射中心对应的真实位置，用d_k表示离真实位置/>

最近的网格点，用δ_k表示真实位置与其最近的网格点之间的距离，则

从而构建出高低频段数据的离网稀疏表征模型：

Y＝ΦX+E

其中，E为噪声矩阵；稀疏矩阵X为高低频段回波相干合成的二维成像结果，X的维度为D×M；离网字典

其维度为N×D；/>

为目标真实位置/>

的所有可能取值；离网字典Φ的每一个原子表示为：

其中，j是复数单位，n为高低频段在合成全频段中的索引向量，其维度为N×1；n₁为低频段在合成全频段中的索引向量；n₂为高频段在合成全频段中的索引向量；

其中d_r为网格点；

稀疏矩阵X中的每一列x都是稀疏向量，其不为零项的值即表示目标散射中心的幅度σ″，不为零项对应的原子相位即对应目标散射中心的位置y_k，于是求解出的稀疏向量x即可表示目标一维距离像，而所有回波求解出的稀疏矩阵X即为高低频段回波相干合成的二维成像结果；

子步骤2.2，由于δ_r取值很小，将原子

在d_r处进行二阶泰勒展开，即

其中

若记基字典矩阵为A＝[a(d₁)，…，a(d_D)]，一阶补偿字典矩阵B＝[b(d₁)，…，b(d_D)]，二阶补偿字典矩阵C＝[c(d₁)，…，c(d_D)]，其中δ＝[δ₁，…，δ_D]^T称为模型补偿向量，于是得到基于二阶泰勒展开的双频段离网稀疏表征模型：

Y＝ΦX+E＝(A+Bdiag(δ)+Cdiag(δ²))X+E

其中

表示两个向量的hadamard积，diag(δ)表示以δ为对角线元素的对角矩阵；根据高低频段在全频段的索引向量，生成对应的基字典矩阵A、一阶补偿字典矩阵B、二阶补偿字典矩阵C；

步骤3，利用稀疏贝叶斯学习算法依次对不同方位单元中所包含的离网散射中心位置和散射强度进行求解，即可实现目标的二维高分辨成像；

步骤3包含以下子步骤：

子步骤3.1，第m次回波的基于二阶泰勒展开的双频段离网稀疏表征模型为：

y＝Φx+ε＝(A+Bdiag(δ)+Cdiag(δ²))x+ε

其中，观测向量y为高低频段观测数据矩阵Y的某一列，稀疏向量x为稀疏矩阵X的某一列；设噪声向量ε的实部和虚部均独立同分布于均值为0、协方差为

的高斯分布，其中β为噪声精度，I为单位矩阵，于是ε就服从于均值为0、协方差为β^-1I的复高斯分布，可以得到已知β的条件下ε的概率密度函数p(ε|β)为：

p(ε|β)＝CN(ε|0，β^-1I)

其中，CN()表示复高斯分布的概率密度函数，设噪声精度β服从于参数为(c，d)的伽马分布，于是得到β的概率密度函数p(β)为：

p(β)＝Gam(β|c，d)

其中，Gam()表示伽马分布的概率密度函数，c、d分别表示伽马分布的形状参数和逆尺度参数；

观测向量y的似然概率p(y|x，β，δ)服从均值为Φx、方差为β^-1I的复高斯分布，即：

p(y|x，β，δ)＝CN(Φx，β^-1I)

将散射中心位置和幅度向量建模为两层的高斯逆伽马先验，即认为x服从于方差为α的复高斯分布，于是得到已知α的条件下x的概率密度函数p(x|α)为：

p(x|α)＝CN(0，Λ^-1I)

其中，方差矩阵Λ＝diag(α)，α＝[α₁，…，α_r，…，α_D]^T∈R^D，并且α_r独立同分布于参数为(a，b)的逆伽马分布，即α_r的概率密度函数p(α_r)为：

p(α_r)＝IGam(α_r|a，b)，(r＝1，...，D)

其中，IGam()表示逆伽马分布的概率密度函数，a、b分别表示逆伽马分布的形状参数和尺度参数；假设补偿向量δ_r服从于均匀分布，即δ_r的概率密度函数p(δ_r)表示为：

子步骤3.2，利用期望最大化EM框架对未知参数进行贝叶斯推断，基于证据最大化理论估计参数θ＝(α，β，δ)，估计值

表示为：

其中，p(θ|y)为已知y的条件下θ的概率密度函数，p(y，θ)为y，θ的联合概率密度函数，即最大化(y，α，β，δ)的联合概率密度函数的对数，其表达式为：

其中，D＝β^-1I+ΦAΦ为边缘似然分布的协方差矩阵；

子步骤3.3，利用期望最大化EM框架通过最大化代价函数的下界Q(α，β，δ)求参数，即

maxQ(α，β，δ)＝maxE_{x～p(x|y，α，β，δ)}{ln p(y，x，α，β，δ)}

下界Q(α，β，δ)等于完全数据的对数联合分布函数ln p(y，x，α，β，δ)对隐变量x的条件概率分布p(x|y，α，β，δ)的期望；根据概率图模型可写出其联合分布为：

p(y，x，α，β，δ)＝p(y|x，β，δ)p(x|α)p(α)p(β)p(δ)

其中，p(δ)为δ的概率密度函数；

E步：求隐变量x的后验期望，x的后验概率分布为p(x|y，α，β，δ)∝p(y|x，β，δ)p(x|α)

＝CN(x|μ，∑)

将期望值μ作为x的估计值，∑是x的后验协方差矩阵；

M步：通过最大化E{log p(x|α)p(α)}得到方差α的更新公式，通过最大化E{log p(y|x，β，δ)p(β)}得到噪声精度β的更新公式，通过最大化E{log p(y|x，β，δ)p(δ)}得到补偿因子δ的更新公式；

子步骤3.4，利用所述期望最大化EM框架依次对不同方位单元中所包含的散射中心位置和散射强度进行求解，即可得到目标的二维ISAR图像。

2.根据权利要求1所述的离网稀疏贝叶斯学习的双频段ISAR成像方法，其特征在于，子步骤3.4包括以下子步骤：

子步骤3.4.1，初始化超参数a＝b＝10^-4，模型补偿向量δ＝{0}_D×1，模型稀疏度为K_m，最大迭代次数Iter_max＝3000，迭代收敛门限τ＝10^-2；

子步骤3.4.2，对于高低频段观测数据矩阵的第m列，初始化噪声精度

其中/>

S_i为子频段i未截取前的二维像，n_o∈I_noise表示噪声所在的距离单元索引，N_i为高低频段距离向单元数，并初始化为c＝1，d＝Nβ_m；初始化稀疏向量x的方差α＝|A^Hy|/N；当前迭代次数iter＝1；

子步骤3.4.3，根据如下公式更新稀疏向量x：

μ＝β∑Φ^Hy

∑＝(βΦ^HΦ+Λ^-1)^-1

将期望值μ作为x的估计值；

子步骤3.4.4，根据如下公式更新α_r：

其中，<|x_r|²>＝E{|x_r|²}＝|μ_r|²+Re(∑_rr)，μ_r为μ的第r个元素，∑_rr为∑的第r行第r列元素，| |表示取模值操作，Re()表示取实部操作；

根据如下公式更新噪声精度β：

其中，

子步骤3.4.5，找到前K_m个最大的α对应的索引I_k，并且通过求方程f′(δ_k)(k∈I_k)的根更新每一个δ_k；其中，f′(δ_k)通过最大化E{log p(y|x，β，δ)p(δ)}推导出来的：

其中

v₁＝Re[diag(μ^*)B^H(y-Aμ)]-Re[diag(B^HA∑)]

v₂＝Re[diag(μ^*)C^H(y-Aμ)]-Re[diag(C^HA∑)]

()^*表示取共轭操作，()^H表示取共轭转置操作，()_k，k表示取出矩阵的第k行第k列元素的操作，

表示取出矩阵的第k行除掉第k列外其他元素的操作，/>

表示取出矩阵的第k列除掉第k行其他元素的操作，()^-k表示取出向量除第k个元素外其余元素的操作，()_k表示取出向量第k个元素的操作；

子步骤3.4.6，判断迭代终止条件

或者iter＜Iter_max是否满足，若有一个条件满足则停止迭代；否则，令迭代次数iter增加1，转到子步骤3.4.3；

子步骤3.4.7，迭代完成后，利用模型量化误差求解目标离网位置

其中/>

表示索引I_k对应的补偿量；然后构造离网字典/>

利用子步骤3.4.3求得的μ作为第m次回波的高分辨一维距离像；

子步骤3.4.8，判断当前方位单元序号是否等于M，若是，输出高分辨二维像；否则，执行子步骤3.4.2。

Images