CN112859075A - 多频带isar融合高分辨成像方法 - Google Patents

Abstract

本发明公开了一种多频带ISAR融合高分辨成像方法，涉及图像处理方法技术领域。所述方法包括如下步骤：首先，直接利用观测数据构造相干字典，建立基于稀疏表示的非相干相位估计模型；采用正交匹配追踪OMP算法求解参数，利用参数对应关系分别估计线性相位和固定相位，并进行非相干相位补偿；建立多频带ISAR融合成像稀疏表示模型，假设散射点系数服从复Gaussian分层先验分布，噪声服从复Gaussian先验分布，直接在复数域利用变分贝叶斯期望最大VB‑EM算法进行求解，得到融合的ISAR高分辨二维图像。所述方法不仅在网格失配情况下能实现较好的非相干相位补偿，而且在低信噪比条件下仍能得到高质量ISAR融合图像。

Description

多频带ISAR融合高分辨成像方法

技术领域

本发明涉及图像处理方法技术领域，尤其涉及一种多频带ISAR融合高分辨成像方法。

背景技术

逆合成孔径雷达(Inverse Synthetic Aperture Radar，ISAR)高分辨图像可以较好地反映目标的二维结构信息，有利于目标的识别与分类，在军事和民用方面应用广泛。多频带ISAR融合成像技术利用在同一视角工作在不同频带的雷达对目标进行观测，得到的回波数据通过信号处理方法融合成一个等效的更大带宽的回波，进而提高目标的距离分辨率，以获得更丰富的目标信息。

当多部雷达同时观测目标时，回波信号往往是不相干的，其相位的非相干量可看作由一个线性相位项和一个固定相位项组成，对多部雷达信号进行相干处理是实现多频带融合成像的前提。目前，非相干相位估计方法主要可分为两大类。一类方法是先将各子频带的频谱进行外推得到重叠频带，再利用最小均方误差准则、最小熵准则或距离像相关等方法估计非相干相位。该类方法原理简单，易于实现，但频带外推误差易受噪声和外推长度的影响，导致非相干相位估计不准，限制了算法的应用范围。另一类方法基于全极点模型，利用求根多重信号分类(Root-Multiple Signal Classification,Root-MUSIC)等谱估计方法估计模型参数，进而得到非相干相位估计值。该类方法虽避免了频带外推，但对全极点模型阶数的估计准确度要求高，且对噪声较为敏感。相比于传统的相干处理方法，基于稀疏表示的相干处理方法不需要重叠频带和估计模型阶数，能够得到更高精度的参数估计，且算法抗噪性能强，估计更加稳健。

经过相干处理后，研究信号融合方法是实现多频带融合高分辨成像的关键。目前，多频带融合成像方法主要可分为基于谱估计和基于稀疏表示两类。谱估计类方法将多频带融合成像问题转化为谱估计问题，利用修正的Root-MUSIC 方法、旋转不变估计信号参数方法(Estimation of Signal Parameter via Rotational Invariance Techniques,ESPRIT)和缺失自适应幅相估计方法(Gapped Adaptive Amplitude and PhaseEstimation,GAPES)等进行参数估计，实现融合成像。该类方法虽然估计精度较高，但需要准确估计散射点个数，这在实际处理中较为困难，且算法性能受噪声水平和数据缺失比的影响较大。稀疏表示类方法利用 ISAR成像信号的稀疏性，将多频带融合成像问题转化为稀疏表示模型，利用稀疏重构方法实现融合成像。有现有技术提出了基于稀疏成分分析的多频带雷达信号综合方法，实现了雷达一维成像，但该方法需要合理选择正则化参数，且运算复杂度较大。还有现有技术从理论上证明了稀疏贝叶斯学习(Sparse Bayesian Learning,SBL)算法实现多频带融合的有效性，利用Gaussian先验进行建模，在求解时无需参数调整且估计精度较高，但模型的灵活性不强。还有现有技术提出了基于Gaussian分层先验的多频带融合融合成像方法，利用两层先验提高了模型的灵活性，但求解过程是在实数域进行，需将复数雷达信号的实部和虚部分别求解，导致运算量较大，而且在实际应用中，复数信号实部和虚部往往具备相同的支撑基或相关性，利用实部和虚部拼接得到的重构结果则无法利用复数信号的这一特性，算法的重构性能有进一步提升的空间。

发明内容

本发明所要解决的技术问题是如何提供一种不仅在网格失配情况下能实现较好的非相干相位补偿，而且在低信噪比条件下仍能得到高质量ISAR融合图像的方法。

为解决上述技术问题，本发明所采取的技术方案是：一种多频带ISAR融合高分辨成像方法，其特征在于包括如下步骤：

首先，直接利用观测数据构造相干字典，建立基于稀疏表示的非相干相位估计模型；

采用正交匹配追踪OMP算法求解参数，利用参数对应关系分别估计线性相位和固定相位，并进行非相干相位补偿；

建立多频带ISAR融合成像稀疏表示模型，假设散射点系数服从复Gaussian 分层先验分布，噪声服从复Gaussian先验分布，直接在复数域利用变分贝叶斯期望最大VB-EM算法进行求解，得到融合的ISAR高分辨二维图像。

采用上述技术方案所产生的有益效果在于：该方法不仅提高了网格失配条件下的相干配准精度，而且在低信噪比和低数据率条件下均能较好地实现多频带高分辨融合成像。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1是本发明实施例中基于复Gaussian分层先验的概率图模型；

图2是本发明实施例所述方法的流程图；

图3a是本发明实施例中非相干距离像结果图；

图3b是本发明实施例中相干补偿后距离像(Q＝1)结果图；

图3c是本发明实施例中相干补偿后距离像(Q＝10)结果图；

图4a是本发明实施例中非相干相位补偿前子频带与全频带频率响应图；

图4b是本发明实施例中非相干相位补偿后子频带与全频带频率响应图；

图5a是本发明实施例中不同信噪比条件下线性相位估计的RMSE；

图5b是本发明实施例中不同信噪比条件下固定相位估计的RMSE；

图6a是本发明实施例中全频带和融合频带频率响应图；

图6b是本发明实施例中子频带、全频带和融合频带一维距离像图；

图7a是本发明实施例中子频带2RD成像图；

图7b是本发明实施例中非相干子频带RD融合成像图；

图8a是本发明实施例中子频带间非相干相位补偿前的频谱图；

图8b是本发明实施例中子频带间非相干相位补偿后的频谱图；

图9a是本发明实施例中RD融合成像图；

图9b是本发明实施例中CGSM融合成像图；

图10a是OMP算法(SNR＝20dB)条件下的融合成像结果图；

图10b是l₁范数稀疏优化算法(SNR＝20dB)条件下的融合成像结果图；

图10c是CGSM算法(SNR＝20dB)条件下的融合成像结果图；

图10d是OMP算法(SNR＝10dB)条件下的融合成像结果图；

图10e是l₁范数稀疏优化算法(SNR＝10dB)条件下的融合成像结果图；

图10f是CGSM算法(SNR＝10dB)条件下的融合成像结果图；

图10g是OMP算法(SNR＝0dB)条件下的融合成像结果图；

图10h是l₁范数稀疏优化算法(SNR＝0dB)条件下的融合成像结果图；

图10i是CGSM算法(SNR＝0dB)条件下的融合成像结果图；

图11a是OMP算法(ρ＝50％)的融合成像结果；

图11b是l₁范数稀疏优化算法(ρ＝50％)的融合成像结果；

图11c是CGSM算法(ρ＝50％)的融合成像结果；

图11d是OMP算法(ρ＝25％)的融合成像结果；

图11e是l₁范数稀疏优化算法(ρ＝25％)的融合成像结果；

图11f是CGSM算法(ρ＝25％)的融合成像结果；

图11g是OMP算法(ρ＝15.6％)的融合成像结果；

图11h是l₁范数稀疏优化算法(ρ＝15.6％)的融合成像结果；

图11i是CGSM算法(ρ＝15.6％)的融合成像结果。

具体实施方式

本发明实施例公开了一种多频带ISAR融合高分辨成像方法，所述方法包括如下步骤：

首先，为避免频带外推和减少网格失配影响，直接利用观测数据构造相干字典，建立基于稀疏表示的非相干相位估计模型，在不增加字典维数的情况下细化网格尺寸，采用正交匹配追踪(Orthogonal Matching Pursuit,OMP)算法求解参数，利用参数对应关系分别估计线性相位和固定相位，并进行非相干相位补偿。

然后，建立融合成像稀疏表示模型，为增强模型灵活性并充分利用噪声的先验信息和散射点分布统计信息，假设散射点系数服从复Gaussian分层 (Complex GaussianScale Mixure,CGSM)先验分布，噪声服从复Gaussian先验分布，为利用复数信号实部与虚部间相同的支撑基特性，直接在复数域利用变分贝叶斯期望最大(Variational BayesianExpectation Maximization,VB-EM) 法进行求解，得到融合的ISAR高分辨二维图像。

进一步的，如图2所示，所述方法主要包括如下步骤：

1)对各子频带回波进行平动补偿等预处理，得到回波的距离频域信号，如下式(1)所示；

2)估计线性相位和固定相位，对各子频带回波进行非相干相位补偿，得到相干的回波信号如下式(6)所示；

3)分别对各子频带进行越分辨单元徙动校正和方位向FFT，得到方位压缩后的回波信号，如下式(9)，得到待融合的观测数据S；

4)设定初始迭代次数g＝1，总迭代次数G＝50，初始化参数a＝b＝c＝d＝10^-4，β₀＝1/var(S)，λ₀＝1/|A₀|，A₀＝Θ^HS，设置收敛门限eps；

5)利用VB-EM法逐脉冲回波进行数据融合，根据如下式(22)、下式(25) 和下式(26)分别更新

直到处理完所有的M个脉冲回波数据，再利用下式(27)全局更新β^(g+1)，即完成一次迭代；

6)判断是否满足终止迭代条件，当

或迭代次数达到设定值G时终止迭代，得到融合后的ISAR二维图像

否则转到步骤5)继续下一次迭代。

下面结合具体内容对上述内容进行说明：

多频带ISAR融合成像回波建模：

假设多部雷达邻近配置，各雷达回波可近似为同一观测视角所得。本申请以工作在不同频带的两个雷达为例进行ISAR融合成像分析，假设雷达均发射 M个脉冲的线性调频(Linear Frequency Modulation，LFM)信号，脉冲重复周期均为T_r，慢时间可表示为t_m＝mT_r(m＝0,1,2,…,M-1)，子频带1和子频带2的载频分别为f_c1和f_c2，工作带宽分别为B₁和B₂。子频带回波经解线频调和平动补偿等预处理后，频域可表示为

其中，P为目标散射点个数，a_p为散射点p的散射系数，c为电磁传播速率，ΔR_p(t_m)为散射点p到参考点的相对距离。对于子频带1，i＝1，f₁为子频带1回波的频率序列，可离散表示为f₁＝f₀+n₁Δf，其中f₀＝f_c1-B₁/2为子频带1的起始频率，Δf为频率采样间隔，n₁＝0,1,…,N₁-1，N₁＝B₁/Δf为子频带1的频率采样点数。对于子频带2，i＝2，f₂为子频带2回波的频率序列，有f₂＝f₀+n₂Δf， n₂＝N-N₂,N-N₂+1,…,N-1，N为全频带的频率采样数，N₂＝B₂/Δf为子频带2的频率采样点数。假设f₀也为全频带的起始频率，则子频带1信号和子频带2信号可看作是从全频带信号中稀疏采样得到的两段信号，有N≥N₁+N₂。多频带融合 ISAR成像即利用观测的子频带1和子频带2的回波数据融合得到全频带数据，通过等效增加信号带宽来提高距离分辨率，实现ISAR高分辨成像。

相干处理

由于发射不同频带信号的雷达之间存在系统时间同步误差、初相差以及目标与各雷达之间的距离差等，各子频带回波信号往往是不相干的。因此，对频带不同的多部雷达回波数据进行相干处理是实现多频带雷达融合成像的前提。

考虑到雷达系统自身硬件差异等原因可能引起的时间延迟和初始相位，则子频带回波的频域可重新写为

其中，τ_i和

分别为子频带中的时延和初始相位。当N₁＝N₂时，子频带1与子频带2的频率采样点数相同，则有n₂＝n₁+N-N₂，f₂＝f₀+n₁Δf+ΔB＝f₁+ΔB，其中，ΔB＝(N-N₂)Δf。子频带2的回波可写为

s₂(n₁,t_m)＝s₁(n₁,t_m)exp(jn₁α+jη) (3)

由式(3)可以看出，不同频带的两个雷达回波之间的非相干相位可看作由一个线性相位项α和一个固定相位项η组成，相干处理的目的即求出这两个非相干相位并实现相位补偿，保证频带回波信号间的相干性，为后续进行多频带融合成像奠定基础。需要指出的是，上述推导均是以两部雷达的频率采样点数相同为前提进行的，即N₁＝N₂。若N₁≠N₂时，可取两者中较小的值N_C＝min(N₁,N₂)，在两段雷达回波中均取出长度为N_C的数据进行相干处理。

非相干相位估计：

为便于求解，只取某一脉冲回波的信号进行分析，由于慢时间 t_m＝mT_r(m＝0,1,2,…,M-1)，对于第m个脉冲回波，式(3)可重写为

其中，

由于线性相位项α的取值范围为[0,2π)，故可将α进行离散化处理，令α＝2πk/K(k＝0,1,…,K-1且K＞＞N₁)。式(4)可稀疏表示为s₂＝Fa，其中，s₂为子频带2第m个脉冲回波信号，可表示为s₂＝[s₂(0,m),s₂(1,m),…,s₂(N₁-1,m)]^T；a为系数矢量，可表示为

F为相干处理字典，可表示为 F＝[F₀,F₁,…F_K-1]，其列原子F_k定义为

当线性相位项α正好位于网格中心时，此时不存在网格失配的情况，利用稀疏表示求解方法可以得到较为准确的估计值。但当线性相位项α不位于网格中心而偏离网格时，利用稀疏表示方法求解时只能估计到与真实值最接近的网格中心上，此时存在网格失配问题，导致线性相位项的估计精度受到影响，进而影响后续估计固定相位项。为了减轻网格失配带来的影响，需要细化网格尺度使得字典更加精细，以提高相位估计精度。若直接增大字典维数K，虽能提高估计精度，但同时会因字典维数过大而增加运算量。为在提高线性相位离散精细化程度的同时不增加字典维数，选择缩小网格区域，即将线性相位α的取值范围由[0,2π)缩小为[0,2π/Q)，其中Q为网格区域缩小参数且Q≥1。此时相干字典中的改进原子

可表示为

(当Q＝1 时，相干字典原子为F_k)。多频带雷达信号相干处理问题可稀疏表示为

由于两部雷达之间的非相干相位是唯一的，即式(5)中向量a中只有一个非零元素，其余元素均为零，满足稀疏性，因此，可利用稀疏重构算法估计非相干相位。所述方法采用OMP算法进行求解，得到系数向量估计值

假设

中模值最大的元素(即非零元素)为第i个元素，对应的值为

则线性相位估计值为

固定相位的估计值为

得到线性相位估计值

和固定相位估计值

后，对子频带2回波进行相位补偿，得到与子频带1回波信号相干的子频带2回波信号

依次对所有的脉冲回波进行上述的非相干相位补偿即可得到相干的频带信号，为后续的多频带ISAR融合成像奠定基础。

基于CGSM的多频带融合ISAR高分辨成像：

基于稀疏表示的融合成像方法比谱估计类融合成像方法具有更好的成像性能和抗噪性能。考虑到雷达回波为复数信号，不同于经典SBL方法中的Gaussian 先验模型，本申请采用复Gaussian分层先验进行建模，增强了模型的灵活性，直接在复数域进行贝叶斯推理，利用VB-EM方法实现求解，得到融合的ISAR 高分辨图像。

融合成像稀疏表示模型：

经过非相干相位补偿后，得到相干的各频带回波，可统一表示为

其中，ΔR_p(t_m)＝y_p cosΔθ_m+x_p sinΔθ_m，(x_p,y_p)为散射点p的坐标，Δθ_m为累积转角，目标可看作匀速转动，假设转动角速度为ω，有Δθ_m＝ωt_m。由于成像过程中累积转角Δθ_m较小，近似有cosΔθ_m≈1，sinΔθ_m≈Δθ_m，则式(6)可写为

经过越分辨单元徙动校正后，子频带回波距离频域可表示为

其中，a′_p＝a_p exp(-j4πf₀ y_p/c)。在利用传统的距离-多普勒(Range-Doppler)算法成像时，一般是先进行距离向逆傅里叶变换得到一维距离像，再进行方位向傅里叶变换得到ISAR二维图像。考虑到多频带融合和距离向压缩过程中存在一定的误差，会对后续的方位向压缩产生影响，故在多频带融合成像时，选择先进行方位向傅里叶变换实现方位压缩，再在距离维进行多频带融合成像，得到二维图像。方位向傅里叶变换变换后得到的待融合子频带信号为

其中，

f_d为多普勒频率。令

由于ω_p∈(0,1]，将数字频率离散化，有ω_p＝l/L，l＝0,1,…,L-1，且有L≥N。

雷达发射M次脉冲后，假设S₁为N₁×M维的子频带1回波数据，S₂为N₂×M维的子频带2回波数据，

为(N₁+N₂)×M维多频带观测回波数据。考虑到实际有噪声的情况，方位向压缩后的多频带融合ISAR二维成像模型可表示为

S＝ΦΨΑ+ε＝ΘΑ+ε (10)

式中，Φ为(N₁+N₂)×N维观测矩阵，Ψ为N×L维字典矩阵，Θ＝ΦΨ为传感矩阵，ε为噪声向量，散射点系数矩阵A即为融合后的目标图像。观测矩阵Φ和字典矩阵Ψ可表示为

分别对每个脉冲回波进行频谱合成，第m次脉冲回波对应的多频带融合成像模型可写为

其中，s₁为该脉冲对应的子频带1回波数据，s₂为该脉冲对应的子频带2 回波数据，Α_·m为该脉冲数据对应的融合结果，且有

稀疏先验模型

为增强模型灵活性，充分利用噪声先验信息和散射点分布统计特性，本申请采用由复Gaussian先验和Gamma先验组合的分层先验进行稀疏先验建模。

假设噪声ε服从均值为0，方差为β^-1的复Gaussian分布，即

则回波S的似然函数也服从复高斯分布，可写为

其中，I为单位矩阵。为方便进行贝叶斯推理，再假设噪声参数β服从与 Gaussian分布共轭的Gamma分布，即

p(β)＝Gamma(β|a,b) (16)

其中，Gamma(β|a,b)＝Γ(a)^-1baα^a-1e^-bβ，

为保证先验的无信息性，a,b一般设置为较小值，如a＝b＝10^-4。

假设散射点系数矩阵A服从分层复高斯先验，即A中各元素A_lm服从均值为 0，方差为λ_lm ^-1的复Gaussian分布，系数之间是相互独立分布的，则散射点系数矩阵A的条件概率密度函数为：

再对超参数λ_·m添加一层相互独立的Gamma分布，则λ_·m的概率密度函数为：

为保证先验的无信息性，c,d一般设置为较小值，如c＝d＝10^-4。为直观地表示稀疏先验建模过程，图1给出了概率图模型，其中，

表示已知的观测回波数据，○表示未知变量，

表示超参数。此时的散射点系数服从由复Gaussian 分布与Gamma分布组成的两层先验分布，这样的分层复Gaussian分布概率模型能比单纯的复Gaussian分布得到更稀疏的解。

基于VB-EM算法的求解

稀疏表示问题求解时，常用的稀疏重构方法主要包括以OMP算法为代表的贪婪类算法、以l₁范数稀疏优化算法为代表的凸优化类算法和以SBL算法为代表的贝叶斯推断类算法。贪婪类算法虽原理简单，但算法精度不高，且易受噪声水平等因素影响，相对来说凸优化类算法精度较高，但易陷入局部最优，且需要人工调整参数。贝叶斯推断类算法能自适应学习参数，利用完全贝叶斯推理避免了人工调整参数，在保证重构精度的同时提高了算法自适应性。因此，本申请采用VB-EM方法进行模型求解。

由于散射点系数矩阵A的每一列A_·m是相互独立的，故可以分别对每个脉冲回波数据S_·m进行重构。变分贝叶斯方法假设各未知变量的后验概率相互独立，则联合后验概率密度可因式分解为

p(A_·m,λ_·m,β|S_·m)≈q(A_·m)q(λ_·m)q(β) (19)

其中，q(·)表示后验概率密度估计。

对于A_·m，其最大后验估计(Maximum a Posterior,MAP)为

其中，<·>表示求期望值，Λ_·m＝diag(<λ_1m>,<λ_2m>,…,<λ_Lm>)表示由超参数λ_lm(l＝1,2,…,L)的期望值组成的对角矩阵；q(A_·m)可看作近似服从均值为μ_·m，协方差为Σ_m的复Gaussian分布，即有q(A_m)～CN(A_m|μ_m,Σ_m)，其中，

Σ_m＝(<β>Θ^HΘ+Λ_m)^-1 (22)

此时所求的均值μ_·m即表示该脉冲回波数据S_·m所对应的目标图像估计值

则融合得到的目标图像为

为得到融合图像

则需要估计尺度参数λ和噪声参数β，。

当先验与似然函数互为共轭时，后验概率密度与先验概率有相同的形式。由式(18)可知尺度参数λ_·m的先验为Gamma分布，与式(17)所表示的似然函数Gaussian分布共轭，则其近似后验概率密度q(λ_·m)也服从Gamma分布，即

其中，

Σ_m-ll表示矩阵Σ_m对角线上第l个元素值，l＝1,2,…,L。

同理，噪声参数β的先验为Gamma先验，与其Gaussian似然函数共轭，则其近似后验概率密度q(β)同样服从Gamma分布，即

其中，

在完全贝叶斯推理中，后验概率的期望一般作为未知变量的估计，即可利用后验概率密度q(A_·m)、q(λ_·m)和q(β)的期望值得到未知变量的估计值。由于q(A_·m) 服从复Gaussian分布，q(λ_·m)和q(β)均服从Gamma分布，可以得到相应的估计值为

<A_·m>＝μ_·m＝<β>Σ_mΘ^HS_·m (25)

利用式(25)、式(26)和式(27)可分别实现对ISAR融合图像A、尺度参数λ和噪声参数β的迭代更新。

仿真实验

利用多频带一维雷达信号和二维ISAR成像实验分别验证所提相干补偿以及融合成像算法的有效性和优越性。

一维信号相干补偿和融合算法性能验证：

雷达系统的参数设置如表1所示。

表1雷达系统参数设置

假设目标距离向上有5个理想散射点，理想的频率响应为

其中，σ_p为散射点p的散射系数，有σ₁＝3,σ₂＝4,σ₃＝5,σ₄＝6,σ₅＝7，f₀＝10GHz为起始频率，Δf＝5MHz为频率采样间隔，ΔR_p为散射点p到参考点的相对距离，有ΔR_p1＝0.5m,ΔR_p2＝0.7m,ΔR_p3＝1m,ΔR_p4＝2m,ΔR_p5＝2.5m。当i＝1时，n₁＝0,1,…,N₁-1，其中 N₁＝100，s₁为子频带1对应的频率响应；当i＝2时，n₂＝N-N₂,N-N₂+1,…,N-1，其中N₂＝100，N＝400，s₂为子频带2对应的频率响应；当n＝0,1,…,N-1时，s为全频带的频率响应。

相干补偿算法性能验证

为验证算法的相干补偿性能，在子频带2的频率响应中添加α＝π/7的线性相位和β＝π/6的固定相位，以子频带1为基准，利用子频带1和子频带2测得的目标频率响应数据进行相干处理。子频带1和子频带2的一维距离像如图3a所示，由于存在非相干相位，其距离像发生了移位。利用所提相干处理法进行非相干相位估计，当K＝300时，若利用原子F_k构造相干处理字典，线性相位可离散为α＝2πk/K，计算得到线性相位真实值对应的列数k＝150/7，并不为整数，故在字典中没有与理想线性相位对应的列原子，只能选择与其最接近的列原子得到线性相位估计值，此时存在网格失配问题。利用改进原子

构造相干处理字典，分别令Q＝1、Q＝10，当Q＝1时对应原子F_k，当Q＝10时对应改进原子

得到的相位估计结果如表2所示。可以看出，Q＝1时估计的线性相位和固定相位相对误差均比Q＝10时大，说明网格失配情况下利用改进原子构造相干处理字典能有效减小线性相位估计误差和固定相位估计误差。为直观地表示相干处理效果，分别利用估计结果对子频带2的频率响应进行相干补偿。Q＝1时子频带1和补偿后子频带2的一维距离像如图3b所示，此时两个频带的一维距离像之间仍存在失配现象，表明此时相位估计有误差，导致非相干量并没有被完全补偿；Q＝10 时子频带1和补偿后子频带2的一维距离像如图3c所示，此时两个子频带的一维距离像基本重合，表明此时两个子频带信号已为相干信号，非相干量已被较好地补偿。这说明利用改进原子构造相干处理字典能有效减小网格失配影响，提高相位估计精度。

表2网格失配情况下相位估计结果

为进一步验证算法相干补偿有效性，利用改进的相干处理字典进行相干补偿。补偿前子频带和全频带的频率响应实部如图4a所示，可以看出由于子频带 2中存在非相干相位，其频率响应与全频带频率响应不匹配。利用所提方法进行相干补偿，补偿后子频带和全频带的频率响应实部如图4b所示，可以看出补偿后的子频带2频率响应与全频带频率响应基本吻合，使得子频带信号相参，验证了相干算法的有效性，为后续进行多频带融合处理奠定了基础。

相干补偿算法抗噪性能验证

为验证算法的抗噪性能，在子频带的频率响应中添加高斯白噪声，SNR步长为3dB，范围为0dB到30dB。分别利用本申请所提相干处理方法与基于 root-music的全极点相干处理方法进行相干补偿，在每个固定的SNR条件下进行100次独立的蒙特卡洛仿真，相位估计的均方根误差(Root Mean Square Error，RMSE)随SNR变化曲线如图5a-5b所示。

从图5a中可以看出，随着SNR减小，利用本申请所提方法估计线性相位项的RMSE基本没有变化，说明在估计线性相位时该算法的鲁棒性较好，而利用极点相干法估计线性相位项的RMSE较大，且随SNR变化较为明显，这是因为全极点相干法对模型阶数的估计精度要求较高，而利用root-music方法等估计模型阶数时易受噪声水平的影响，导致SNR较低时可能无法正确估计出极点个数，使得线性相位项估计误差较大。从图5b中可以看出，随着SNR减小，两种方法固定相位项的RMSE都有所增大，其中极点相干法的变化幅度较大，特别是SNR低于15dB时估计误差急剧增大，无法正确估计固定相位，算法已基本失效。相比之下，利用本申请所提方法估计固定相位的RMSE随SNR的变化幅度较小，即使在SNR为0dB时仍能较为准确地估计出线性相位和固定相位，体现了所提算法具有更强的抗噪性。

融合算法有效性验证:

在子频带相干补偿后，利用所提融合算法进行子频带融合，全频带和融合频带的频率响应如图6a所示，可以看出融合后的频率响应与全频带基本重合，说明所提算法能有效实现频带融合。利用融合算法得到的融合频带一维距离像以及子频带1和全频带直接FFT得到的一维距离像如图6b所示。通过比较可以看出，子频带1由于带宽较小，分辨率较低，无法区分出5个散射点；全频带由于频带增加，分辨率有所提高，可以基本区分出5个散射点，但由于直接利用FFT得到的一维距离像，受主瓣展宽、能量泄漏等影响，导致ΔR_p1＝0.5m和ΔR_p2＝0.7m对应的两个散射点无法被完全区分；利用所提融合算法得到的一维距离像可以完全区分出5个散射点，这是因为该方法利用稀疏表示理论，突破了传统傅里叶理论方法中带宽Rayleigh限的限制，达到了超分辨的效果。

二维实测数据实验:

实测数据采用Yak-42飞机回波数据进行实验，信号中心频率为5.52GHz，带宽为400MHz，共发射256个脉冲，脉冲重复频率为100Hz，脉冲宽度为25.6μs，采样频率为10MHz，全频带回波数据中共有256个频率采样点。

算法有效性验证:

分别选取全频带回波的前64个和后64个频点作为子频带1和子频带2的数据。为验证相干处理算法的有效性，在子频带2的回波数据中再添加将π/9的线性相位项和π/8的固定相位项，为回波添加Gaussian白噪声，使SNR为20dB。此时子频带2RD成像结果和非相干子频带直接RD融合成像结果分别如图7a、图7b所示，可以看出，由于存在非相干相位，无法正确地实现多频带融合成像。利用本申请所提算法进行相干处理，分别对每个回波脉冲进行非相干相位估计并补偿，图8a-8b给出了第128个脉冲对应的子频带间非相干相位补偿前后的频谱图。从图8a可以看出，由于子频带2中存在非相干相位，两子频带信号之间是非相参的，从图8b可以看出，利用所提算法对子频带2进行相干补偿后，两子频带信号之间是相参的，验证了相干处理算法的有效性。

图9a-图9b给出了相干处理后多频带RD直接融合成像结果和利用本申请算法融合成像结果。对比图9a和图9b可以看出，直接采用RD算法融合成像时由于频带数据缺失导致成像结果中产生能量泄露且无法得到清晰的目标图像结果，而采用本申请所提的CGSM算法进行融合成像时得到清晰的目标图像，验证了所提融合算法的有效性。

融合算法优越性验证:

A不同SNR条件下算法的抗噪性

改变噪声水平，使SNR分别为0dB，10dB，20dB，分别利用OMP算法、 l₁范数稀疏优化算法和本申请所提算法在不同SNR条件下进行多频带融合 ISAR成像，均设置相同的收敛条件，成像结果如图10所示。从图10a、图10d 和图10g中可以看出，随着SNR降低，采用OMP算法得到的融合ISAR图像质量越来越差，特别是当SNR低于10dB时，成像结果中存在明显的虚假散射点，严重影响了图像质量。这是因为OMP算法性能易受噪声水平的影响，抗噪性能不强，在低信噪比条件下对噪声的抑制效果不佳。从图10b、图10e和图10h可以看出，随着SNR降低，采用l₁范数稀疏优化算法得到的融合ISAR 图像质量也有所下降，特别是当SNR为0dB时，融合结果中也存在一定的虚假散射点，但与同SNR条件下采用OMP算法得到的融合结果相比，图像质量有所改善。这是因为l₁范数稀疏优化算法的抗噪性能和重构精度比OMP算法强，但其正则化参数与噪声水平有关，若估计不准也会对算法性能造成影响。从图10c、图10f和图10i可以看出，随着SNR降低，采用本申请所提算法均能得到质量较好的融合ISAR图像，说明本申请算法的抗噪性能较好，在低信噪比条件下仍能较好地抑制噪声。

为更直观地对比融合成像结果的图像质量，采用图像熵作为衡量指标，熵值越小表示图像越聚焦，融合成像结果的图像熵如表3所示。可以看出，在相同SNR条件下，本申请所提方法得到的成像结果熵值最小，其次是l₁范数稀疏优化算法，而OMP算法的成像结果熵值最大，说明本申请所提算法的融合成像质量比其他两种算法都好。另外，随着SNR降低，本申请所提方法的融合图像熵值变化量比其他两种方法都小，说明本申请方法融合成像的抗噪性能比其他两种方法强。

表3不同SNR条件下融合成像结果图像熵

B不同带宽条件下算法的成像性能：

在SNR为10dB时，改变有效数据率ρ，验证不同带宽条件下多频带融合成像性能。在N₁＝N₂＝64(即ρ＝50％)、N₁＝N₂＝32(即ρ＝25％)和N₁＝N₂＝20 (即ρ＝15.6％)三种情况下，分别利用OMP算法、l₁范数稀疏优化算法和本申请所提CGSM算法进行融合成像，结果如图11a-图11i所示。从图11a-11c可以看出，当ρ＝50％时，三种算法均能较好地实现多频带融合成像，得到目标基本轮廓，但利用OMP算法和l₁范数稀疏优化算法得到的成像结果背景不够干净，出现少许的虚假散射点。从图11d-11i可以看出，随着有效数据率ρ降低，利用OMP算法和l₁范数稀疏优化算法的成像性能急剧下降，成像结果中出现较多的虚假散射点，导致无法分辨出目标的几何结构，而利用本申请算法仍能得到较为清晰的目标图像，当ρ＝15.6％时也能较好地重建出目标基本轮廓，说明本申请算法具有较强的稳健性，仅利用较少的有效数据也能实现多频带融合成像。

不同数据率条件下融合成像结果的图像熵如表4所示，相同数据率条件下利用本申请所提方法得到图像的熵值比l₁范数稀疏优化算法得到图像的熵值小，同样说明了利用本申请算法能够得到较好的融合成像结果。

表4不同数据率融合成像结果图像熵

从以上分析可以看出，该方法不仅提高了网格失配条件下的相干配准精度，而且在低信噪比和低数据率条件下均能较好地实现多频带高分辨融合成像。

Claims (6)

1.一种多频带ISAR融合高分辨成像方法，其特征在于包括如下步骤：

首先，直接利用观测数据构造相干字典，建立基于稀疏表示的非相干相位估计模型；

采用正交匹配追踪OMP算法求解参数，利用参数对应关系分别估计线性相位和固定相位，并进行非相干相位补偿；

建立多频带ISAR融合成像稀疏表示模型，假设散射点系数服从复Gaussian分层先验分布，噪声服从复Gaussian先验分布，直接在复数域利用变分贝叶斯期望最大VB-EM算法进行求解，得到融合的ISAR高分辨二维图像。

2.如权利要求1所述的多频带ISAR融合高分辨成像方法，其特征在于，所述方法包括如下步骤：

1)对各子频带回波信号进行预处理，得到回波的距离频域信号；

2)估计线性相位和固定相位，对各子频带回波信号进行非相干相位补偿，得到相干的回波信号；

3)分别对各子频带进行越分辨单元徙动校正和方位向FFT变换，得到方位压缩后的回波信号以及待融合的观测数据S；

4)设定初始迭代次数g，总迭代次数G，初始化参数a＝b＝c＝d＝10^-4，β₀＝1/var(S)，λ₀＝1/|A₀|，A₀＝Θ^HS，并设置收敛门限eps；

5)利用VB-EM法逐脉冲回波进行数据融合，并分别更新协方差

ISAR融合图像

尺度参数

直到处理完所有的M个脉冲回波数据，再全局更新噪声参数β^(g+1)，即完成一次迭代；

6)判断是否满足终止迭代条件，当

或迭代次数达到设定值G时终止迭代，得到融合后的ISAR二维图像

否则转到步骤5)继续下一次迭代。

3.如权利要求2所述的多频带ISAR融合高分辨成像方法，其特征在于，得到的回波的距离频域信号如式(1)所示：

脉冲重复周期为T_r，慢时间可表示为t_m＝mT_r(m＝0,1,2,…,M-1)，子频带1和子频带2的载频分别为f_c1和f_c2，工作带宽分别为B₁和B₂；子频带回波预处理后，频域可表示为：

其中，P为目标散射点个数，a_p为散射点p的散射系数，c为电磁传播速率，ΔR_p(t_m)为散射点p到参考点的相对距离；对于子频带1，i＝1，f₁为子频带1回波的频率序列，可离散表示为f₁＝f₀+n₁Δf，其中f₀＝f_c1-B₁/2为子频带1的起始频率，Δf为频率采样间隔，n₁＝0,1,…,N₁-1，N₁＝B₁/Δf为子频带1的频率采样点数；对于子频带2，i＝2，f₂为子频带2回波的频率序列，有f₂＝f₀+n₂Δf，n₂＝N-N₂,N-N₂+1,…,N-1，N为全频带的频率采样数，N₂＝B₂/Δf为子频带2的频率采样点数；假设f₀也为全频带的起始频率，则子频带1信号和子频带2信号可看作是从全频带信号中稀疏采样得到的两段信号，有N≥N₁+N₂。

4.如权利要求3所述的多频带ISAR融合高分辨成像方法，其特征在于：得到相干的回波信号如式(6)所示：

考虑到时间延迟和初始相位，则子频带回波的频域可重新写为

其中，τ_i和

分别为子频带中的时延和初始相位；当N₁＝N₂时，子频带1与子频带2的频率采样点数相同，则有n₂＝n₁+N-N₂，f₂＝f₀+n₁Δf+ΔB＝f₁+ΔB，其中，ΔB＝(N-N₂)Δf；子频带2的回波可写为

s₂(n₁,t_m)＝s₁(n₁,t_m)exp(jn₁α+jη) (3)

其中：α为线性相位项和η为一个固定相位项；

由于慢时间t_m＝mT_r(m＝0,1,2,…,M-1)，对于第m个脉冲回波，式(3)可重写为

其中，

由于线性相位项α的取值范围为[0,2π)，故可将α进行离散化处理，令α＝2πk/K(k＝0,1,…,K-1且K＞＞N₁)；式(4)可稀疏表示为s₂＝Fa，其中，s₂为子频带2第m个脉冲回波信号，可表示为s₂＝[s₂(0,m),s₂(1,m),…,s₂(N₁-1,m)]^T；a为系数矢量，可表示为

F为相干处理字典，可表示为F＝[F₀,F₁,…F_K-1]，其列原子F_k定义为

将线性相位α的取值范围由[0,2π)缩小为[0,2π/Q)，其中Q为网格区域缩小参数且Q≥1；此时相干字典中的改进原子

可表示为

(当Q＝1时，相干字典原子为F_k)；多频带雷达信号相干处理问题可稀疏表示为：

采用OMP算法进行求解，得到系数向量估计值

假设

中模值最大的元素为第i个元素，对应的值为

则线性相位估计值为

固定相位的估计值为

得到线性相位估计值

和固定相位估计值

后，对子频带2回波进行相位补偿，得到与子频带1回波信号相干的子频带2 回波信号

经过非相干相位补偿后，得到相干的各频带回波，可统一表示为

其中，ΔR_p(t_m)＝y_pcosΔθ_m+x_psinΔθ_m，(x_p,y_p)为散射点p的坐标，Δθ_m为累积转角，目标可看作匀速转动，假设转动角速度为ω，有Δθ_m＝ωt_m。

5.如权利要求4所述的多频带ISAR融合高分辨成像方法，其特征在于，得到方位压缩后的回波信号如式(9)所示，多频带融合ISAR二维成像的稀疏表示模型如式(10)所示：

由于成像过程中累积转角Δθ_m较小，近似有cosΔθ_m≈1，sinΔθ_m≈Δθ_m，则式(6)可写为

经过越分辨单元徙动校正后，子频带回波距离频域可表示为

其中，a′_p＝a_pexp(-j4πf₀y_p/c)；

方位向傅里叶变换变换后得到的待融合子频带信号为：

其中，

f_d为多普勒频率；令

由于ω_p∈(0,1]，将数字频率离散化，有ω_p＝l/L，l＝0,1,…,L-1，且有L≥N。雷达发射M次脉冲后，假设S₁为N₁×M维的子频带1回波数据，S₂为N₂×M维的子频带2回波数据，

为(N₁+N₂)×M维多频带观测回波数据。考虑到实际有噪声的情况，方位向压缩后的多频带融合ISAR二维成像模型可表示为

S＝ΦΨΑ+ε＝ΘΑ+ε (10)

式中，Φ为(N₁+N₂)×N维观测矩阵，Ψ为N×L维字典矩阵，Θ＝ΦΨ为传感矩阵，ε为噪声向量，散射点系数矩阵A即为融合后的目标图像。观测矩阵Φ和字典矩阵Ψ可表示为

6.如权利要求5所述的多频带ISAR融合高分辨成像方法，其特征在于，利用VB-EM法逐脉冲回波进行数据融合，根据式(22)、式(25)和式(26)分别更新

直到处理完所有的M个脉冲回波数据，再利用式(27)全局更新β^(g+1)：

分别对每个脉冲回波进行频谱合成，第m次脉冲回波对应的多频带融合成像模型可写为

其中，s₁为该脉冲对应的子频带1回波数据，s₂为该脉冲对应的子频带2回波数据，Α_·m为该脉冲数据对应的融合结果，且有

假设噪声ε服从均值为0，方差为β^-1的复Gaussian分布，即

则回波S的似然函数也服从复高斯分布，可写为

其中，I为单位矩阵。为方便进行贝叶斯推理，再假设噪声参数β服从与Gaussian分布共轭的Gamma分布，即

p(β)＝Gamma(β|a,b) (16)

其中，Gamma(β|a,b)＝Γ(a)^-1b^aα^a-1e^-bβ，

为保证先验的无信息性，a,b一般设置为较小值，如a＝b＝10^-4。

假设散射点系数矩阵A服从分层复高斯先验，即A中各元素A_lm服从均值为0，方差为λ_lm ^-1的复Gaussian分布，系数之间是相互独立分布的，则散射点系数矩阵A的条件概率密度函数为：

再对超参数λ_·m添加一层相互独立的Gamma分布，则λ_·m的概率密度函数为：

由于散射点系数矩阵A的每一列A_·m是相互独立的，故可以分别对每个脉冲回波数据S_·m进行重构；变分贝叶斯方法假设各未知变量的后验概率相互独立，则联合后验概率密度可因式分解为：

p(A_·m,λ_·m,β|S_·m)≈q(A_·m)q(λ_m)q(β) (19)

其中，q(·)表示后验概率密度估计；

对于A_m，其最大后验估计MAP为

其中，<·>表示求期望值，Λ_·m＝diag(<λ_1m>,<λ_2m>,…,<λ_Lm>)表示由超参数λ_lm(l＝1,2,…,L)的期望值组成的对角矩阵；q(A_·m)可看作近似服从均值为μ_·m，协方差为Σ_m的复Gaussian分布，即有q(A_m)～CN(A_m|μ_m,Σ_m)，其中，

Σ_m＝(<β>Θ^HΘ+Λ_m)^-1 (22)

此时所求的均值μ_·m即表示该脉冲回波数据S_·m所对应的目标图像估计值

则融合得到的目标图像为

为得到融合图像

则需要估计尺度参数λ和噪声参数β；

由式(18)可知尺度参数λ_·m的先验为Gamma分布，与式(17)所表示的似然函数Gaussian分布共轭，则其近似后验概率密度q(λ_·m)也服从Gamma分布，即

其中，

Σ_m-ll表示矩阵Σ_m对角线上第l个元素值，l＝1,2,…,L。

同理，噪声参数β的先验为Gamma先验，与其Gaussian似然函数共轭，则其近似后验概率密度q(β)同样服从Gamma分布，即

其中，

在完全贝叶斯推理中，后验概率的期望一般作为未知变量的估计，即可利用后验概率密度q(A_·m)、q(λ_·m)和q(β)的期望值得到未知变量的估计值；由于q(A_·m)服从复Gaussian分布，q(λ_·m)和q(β)均服从Gamma分布，可以得到相应的估计值为

<A_·m>＝μ_·m＝<β>Σ_mΘ^HS_·m (25)

利用式(25)、式(26)和式(27)分别实现对ISAR融合图像A、尺度参数λ和噪声参数β的迭代更新。

Images