CN112859074B - 多频带多视角isar融合成像方法 - Google Patents

多频带多视角isar融合成像方法 Download PDF

Info

Publication number
CN112859074B
CN112859074B CN202110049776.2A CN202110049776A CN112859074B CN 112859074 B CN112859074 B CN 112859074B CN 202110049776 A CN202110049776 A CN 202110049776A CN 112859074 B CN112859074 B CN 112859074B
Authority
CN
China
Prior art keywords
radar
echo
imaging
isar
matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110049776.2A
Other languages
English (en)
Other versions
CN112859074A (zh
Inventor
刘利民
郭宝锋
朱晓秀
胡文华
史林
马俊涛
薛东方
曾慧燕
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Army Engineering University of PLA
Original Assignee
Army Engineering University of PLA
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Army Engineering University of PLA filed Critical Army Engineering University of PLA
Priority to CN202110049776.2A priority Critical patent/CN112859074B/zh
Publication of CN112859074A publication Critical patent/CN112859074A/zh
Application granted granted Critical
Publication of CN112859074B publication Critical patent/CN112859074B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/904—SAR modes
    • G01S13/9064—Inverse SAR [ISAR]

Abstract

本发明公开了一种多频带多视角ISAR融合成像方法，涉及图像处理方法技术领域。所述方法首先，建立多频带多视角ISAR融合成像的矢量化稀疏表示模型，为了减少求解时的运算复杂度，提高成像效率，在Bregman迭代算法的基础上，结合加权残差回代和基矩阵条件数优化进一步提高算法的迭代收敛速度。仿真实验验证了算法在实现一维稀疏信号重构以及多频带多视角ISAR融合成像时有效性和优越性。

Description

多频带多视角ISAR融合成像方法
技术领域
本发明涉及图像处理方法技术领域，尤其涉及一种多频带多视角ISAR融 合成像方法。
背景技术
逆合成孔径雷达(Inverse Synthetic Aperture Radar，ISAR)可在远距离、 全天候条件下对目标进行高分辨成像，在军事和民用领域都得到了广泛的应用。 对于单ISAR成像系统而言，通常通过提高发射信号带宽和增大观测相干积累 时间来分别提高成像的距离分辨率和方位分辨率，但同时会带来雷达系统硬件 复杂度高、制造成本大、运动补偿困难等问题，导致直接提升的单雷达系统成 像分辨率有限。多频带多视角ISAR融合成像技术利用工作在不同频带从不同 角度对目标观测得到的多部雷达回波，在信号级进行融合，得到一个更大带宽 和更大视角的信号回波，打破了传统单雷达成像分辨率受发射信号带宽和观测 累积转角的约束，可同时提高雷达成像二维分辨率。
谱估计类方法是一类传统的多雷达融合成像方法。此类方法包括状态空间 法、Root-MUSIC算法以及基于旋转不变估计信号参数方法(Estimation of Signal Parametervia Rotational Invariance Techniques,ESPRIT)算法的融合算法等。虽 然此类方法在精确估计散射点个数的前提下参数估计精度高，但在实际情况下 散射点个数一般很难被精确估计，影响了算法性能。另外，在多频带多视角ISAR 二维融合过程中，散射点的二维坐标匹配也是一个较大的难点。
基于稀疏表示的方法是一类新兴的多雷达融合成像方法。它将多频带多视 角ISAR二维融合成像问题归结为一个信号的稀疏表示问题进行求解，主要包 括以正交匹配追踪算法(Orthogonal Matching Pursuit,OMP)为代表的贪婪类重 构算法、以基追踪算法(Basis Pursuit,BP)为代表的凸优化类重构算法以及稀 疏贝叶斯学习算法(SparseBayesian Learning,SBL)为代表的贝叶斯推理类算 法等。此类方法通过将目标场景离散化，在二维网格中模拟可能存在的散射点， 求解时无需坐标配对也无需估计散射点个数，因此参数估计性能要优于现代谱 估计类方法。但在对二维融合成像模型进行稀疏表示时，需要利用行列堆叠的 一维矢量形式进行求解，数据维数较大。此时，若采用贪婪类算法求解，如OMP 算法，重构精度有限，且需要预知信号的稀疏度，这在实际过程中难以准确估计，制约了算法的应用性；若采用贝叶斯推理类算法求解，估计精度较高，但 涉及到大量的求逆步骤，导致算法的运算量很大，不利于实时处理。因此，寻 找有效、简单、快速的重构方法对实现多频带多视角融合成像来说十分重要。
发明内容
本发明所要解决的技术问题是如何提供一种在求解时能够避免大量的矩阵 求逆过程，大大提高成像效率的多频带多视角ISAR融合成像方法。
为解决上述技术问题，本发明所采取的技术方案是：一种多频带多视角 ISAR融合成像方法，其特征在于包括如下步骤：
1)对各雷达回波信号进行预处理以及互相干处理，得到相干的距离频域- 方位慢时间域回波信号；
2)构造字典矩阵，将回波信号离散化表示；
3)分别将回波数据矢量化拼接，得到待融合的观测信号
Figure BDA0002898636170000021
以及对应的基矩 阵
Figure BDA0002898636170000022
4)利用FLBI算法进行迭代求解，得到目标图像矢量估计值
Figure BDA0002898636170000023
5)将一维目标图像矢量估计值
Figure BDA0002898636170000024
转换为二维矩阵
Figure BDA0002898636170000025
即为多视角多频带 雷达信号融合得到的目标图像。
采用上述技术方案所产生的有益效果在于：所述方法在建立多频带多视角 融合成像稀疏表示模型的基础上，直接在复数域利用LBI算法进行求解，同时 利用加权残差回代和优化基矩阵条件数相结合的方式来进一步加快迭代收敛速 度，在求解时避免了大量的矩阵求逆过程，大大提高了成像效率，仿真实验验 证了算法的有效性和优越性。
附图说明
下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1是本发明实施例中多频带多视角双雷达系统二维观测数据融合示意 图；
图2是本发明实施例中基于矢量化处理的多视角多频带双雷达观测数据融 合示意图；
图3是本发明实施例中所述方法的流程图；
图4a是本发明实施例中目标模型；
图4b是本发明实施例中RD成像结果图；
图5a是本发明实施例中雷达1RD成像结果图；
图5b是本发明实施例中雷达2RD成像结果图；
图5c是本发明实施例中RD融合成像结果图；
图5d是本发明实施例中OMP算法融合成像结果图；
图5e是本发明实施例中所述方法融合成像结果图；
图6a是本发明实施例中波音727飞机回波数据一维距离像；
图6b是本发明实施例中波音727飞机回波数据RD成像结果图；
图7a是本发明实施例中波音727飞机的雷达1RD成像结果图；
图7b是本发明实施例中波音727飞机的雷达2RD成像结果图；
图7c是本发明实施例中波音727飞机的RD融合成像结果图；
图7d是本发明实施例中波音727飞机的OMP算法融合成像结果图；
图7e是本发明实施例中波音727飞机通过本申请所述方法融合成像结果 图；
图8a是本发明实施例中M1＝M2＝16，N1＝N2＝64时观测二维回波数据图(情形1)；
图8b是本发明实施例中M1＝M2＝16，N1＝N2＝64时OMP算法融合成像结果图(情 形1)；
图8c是本发明实施例中M1＝M2＝16，N1＝N2＝64时本申请所述方法的融合成像 结果图(情形1)；
图9a是本发明实施例中M1＝M2＝16，N1＝N2＝32时观测二维回波数据图(情形2)；
图9b是本发明实施例中M1＝M2＝16，N1＝N2＝32时OMP算法融合成像结果图(情 形2)；
图9c是本发明实施例中M1＝M2＝16，N1＝N2＝32时本申请所述方法的融合成像 结果图(情形2)；
图10a是本发明实施例中M1＝M2＝8，N1＝N2＝64时观测二维回波数据图(情形 3)；
图10b是本发明实施例中M1＝M2＝8，N1＝N2＝64时OMP算法融合成像结果图(情 形3)；
图10c是本发明实施例中M1＝M2＝8，N1＝N2＝64时本申请所述方法的融合成像 结果图(情形3)。
具体实施方式
下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、 完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全 部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性 劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明 还可以采用其他不同于在此描述的其它方式来实施，本领域技术人员可以在不 违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例 的限制。
总体的，本发明实施例公开了一种多频带多视角ISAR融合成像方法，所 述方法包括如下步骤：
首先，将ISAR回波信号构建为二维稀疏表示模型，在此基础上，利用矢 量化操作将多频带多视角ISAR融合成像问题转换成一维信号矢量的稀疏重构 问题。其次，为避免贝叶斯推理类算法重构时涉及大量的矩阵求逆操作带来的 运算复杂度较大的问题，提出了利用快速线性Bregman迭代算法进行稀疏重构， 结合加权残差回代和基矩阵条件数优化的思想，进一步提高了算法的迭代收敛 速度。仿真实验表明算法能有效提高迭代收敛速度，较好地实现多频带多视角 ISAR融合成像。
进一步的，如图3所示，本发明所述方法包括如下步骤：
1)对各雷达回波进行运动补偿等预处理以及互相干处理，得到相干的距离 频域-方位慢时间域回波信号；
2)构造字典矩阵，将回波信号离散化表示；
3)分别将回波数据矢量化拼接，得到待融合的观测信号
Figure BDA0002898636170000051
以及对应的基矩 阵
Figure BDA0002898636170000052
4)利用FLBI算法进行迭代求解，得到目标图像矢量估计值
Figure BDA0002898636170000053
5)将一维矢量
Figure BDA0002898636170000054
转换为二维矩阵
Figure BDA0002898636170000055
即为多视角多频带雷达信号融合得 到的目标图像。
下面结合具体步骤对以上方法进行详细说明：
多频带多视角ISAR融合成像稀疏表示模型：
在单部雷达系统观测条件下，由于有限的发射信号带宽和观测角度，导致 其距离分辨率和方位分辨率受到限制。而在多雷达系统观测条件下，若能利用 不同工作频带和不同观测视角的多个雷达同时对目标进行观测，再利用信号处 理技术将观测数据融合成一个更大带宽和更大视角的回波信号，则可同时提高 雷达成像的二维分辨率。值得注意的是，为实现融合成像，多部雷达观测目标 同一散射中心的信息不能差距太大，因此雷达间的观测视角差距不能太大。在 高频段，频带较窄、观测角度较小的情况下，散射点的复幅度可认为是常数， 即不同雷达观测到同一目标的散射响应特性可近似为理想散射点模型。假设所 讨论的多雷达信号已经过预处理和非相干量补偿，实现了信号的完全相关，本 申请以两个不同频带不同视角的雷达信号为例进行ISAR融合成像分析。
单雷达系统ISAR成像模型：
假设ISAR发射线性调频信号，在远场条件下，经解线频调处理后，雷达 回波可表示为
Figure BDA0002898636170000061
其中，fm为距离频率，tn＝nTr为慢时间函数，n＝0,1,…N-1，N为回波脉冲数， Tr为脉冲重复时间，P为目标散射点个数，ap为散射点p的散射系数，c为电磁 传播速率，ΔRp(tn)＝yp cos(Δθ(tn))+xp sin(Δθ(tn))为散射点p到参考点的距离，(xp,yp)为 散射点p的坐标，Δθ(tn)为观测累积转角。
由于成像观测时间短，累积转角Δθ(tn)较小，有cos(Δθ(tn))≈1，sin(Δθ(tn))≈Δθn。 经运动补偿后，目标运动模型可近似为转台模型，假设匀速转动的角速度为ω，则有Δθ(tn)＝ωtn，则式(1)可近似为
Figure BDA0002898636170000062
对频率进行离散化采样，令fm＝f0+mΔf，m＝0,1,…M-1，其中，M为频率采样 点数，Δf为频率采样间隔。在有限带宽和小角度的情况下，忽略散射点越分辨 单元徙动(Migration Through Resolution Cells,MTRC)的影响，雷达回波距离频 域可离散表示为
Figure BDA0002898636170000063
其中，a′p＝apexp(-j4πf0yp/c)。若存在MTRC，可利先利用Keystone变换法重 新定义一个慢时间，以校正由快时间频域-慢时间耦合项引起的越距离单元徙动 校正，再构造一个相位补偿项以校正越多普勒单元徙动。
令
Figure BDA0002898636170000064
一般ISAR成像目标的尺寸较小，有ωm,ωn∈(0,1]，将其离散化，令ωm＝k/K(k＝0,1,…,K-1)，ωn＝l/L(l＝0,1,…,L-1)，且有K≥M，L≥N。
此时，式(3)可表示为
Figure BDA0002898636170000065
其中，akl表示每个像素网格的散射系数幅度。式(4)等价于将目标成像场景 离散化，距离向和方位向分别划分为K个和L个网格，一共有K×L个成像像素 网格。当某一网格交点的坐标(l,k)上存在等效散射点时，此网格点的幅度akl≠0， 反之，当此位置上不存在等效散射点时，则akl＝0。由于目标ISAR图像由有限 个散射点组成，仅占成像场景中很小的一部分，即akl中仅有少数幅度为非零， 大多数幅度为零，故ISAR图像满足很强的稀疏性，适用于稀疏表示理论。
构造距离向傅里叶变换矩阵FR＝[FR(0) FR(1) … FR(m) … FR(M-1)]T，其中 FR(m)＝[exp(-j2πm·0/K) exp(-j2πm·1/K) … exp(-j2πm·(K-1)/K)]，构造方位向傅里叶变 换矩阵FA＝[FA(0) FA(1) … FA(n) … FA(N-1)]，其中 FA(n)＝[exp(-j2πn·0/L)exp(-j2πn·1/L) … exp(-j2πn·(L-1)/L)]T，则ISAR回波可矩阵表 示为
S＝FRΑFA (5)
其中，S为大小为M×N的雷达回波数据，FR为大小为M×K的距离向稀疏基 矩阵，FA为大小为L×N的方位向稀疏基矩阵，Α为大小为K×L的散射系数矩阵， 该矩阵可代表目标二维ISAR图像，其中第k行l列元素为akl。
多频带多视角ISAR融合原理：
本申请以两部雷达为例进行多频带多视角ISAR融合成像研究。两部雷达 相近放置，分别工作在不同频带并以不同视角同时对目标进行观测。雷达1发 射信号频带为
Figure BDA0002898636170000071
共有M1个频率采样点，对应观测角度为
Figure BDA0002898636170000072
共有N1个角度；雷达2发射信号频带为
Figure BDA0002898636170000073
共有M2个频率采样点， 观测角度为
Figure BDA0002898636170000074
共有N2个角度。假设Δf和Δθ分别表示频率采样间隔和 角度采样间隔，M和N分别表示全频带全角度的频率采样个数和角度采样个数， 则全频带的频率采样数据可表示为fm＝f0+mΔf(m＝0,1,…,M-1)，全角度的角度采样 数据可表示为θn＝θ0+nΔθ(n＝0,1,…,N-1)。多视角多频带双雷达系统二维观测数据融 合原理如图1所示，其中，左上点代表雷达1观测数据，右下点代表雷达2观 测数据，空白部分为缺失的频带和视角数据。多频带多视角数据融合即利用雷 达1和雷达2观测得到的已知回波数据，通过稀疏信号重构方法重建出全频带 全视角回波数据，补全频带和视角缺失数据，得到一个等效的更宽频带更大视 角的二维雷达回波，从而提高图像的距离向和方位向分辨率。
由于回波信号行列之间具有耦合性，不能简单地通过分别左乘和右乘行、 列采样矩阵来构造二维融合成像的线性系统模型。因此，采用一维矢量化处理 方式建立多频带多视角ISAR融合成像的稀疏表示模型。将回波二维矩阵按照 频率-角度排为一维向量，即将回波S按列矢量化堆叠，可将式(5)的二维信 号模型转化为一维矢量形式
Figure BDA0002898636170000081
其中，s＝vect(S)，a＝vect(A)，vect(·)表示将矩阵各列向量依次堆叠为一个矢量的操作，
Figure BDA0002898636170000082
表示矩阵(FA)T和FR的Kronecker乘积，即F为全频带全 视角回波矢量对应的基矩阵。基于矢量化处理的多频带多视角双雷达观测数据 融合示意图如图2所示。其中红色矩形部分表示与雷达1的N1个观测角度对应 的基，蓝色矩形部分表示与雷达2的N2个观测角度对应的基。
在全频带全视角回波矢量数据s中去除数据缺失的部分即可得到雷达1和 雷达2的观测回波矢量数据
Figure BDA0002898636170000089
同时在基矩阵F中也去除缺失数据对应的行，令 FA1＝[FA(0) FA(1)… FA(N1-1)]，FA2＝[FA(N-N2) FA(N-N2+1) … FA(N-1)]， FR1＝[FR(0) FR(1) … FR(M1-1)]T，FR2＝[FR(M-M2) FR(M-M2+1) … FR(M-1)]T，则
Figure BDA0002898636170000083
表示雷达1观测回波矢量数据对应的基矩阵，维数为M1N1×KL，对 应图2中去除数据缺失部分后的上侧矩形块部分，
Figure BDA0002898636170000084
表示雷达2观 测回波矢量数据对应的基矩阵，维数为M2N2×KL，对应图2中去除数据缺失部 分后的下侧矩形块部分，则
Figure BDA0002898636170000085
表示双雷达观测回波矢量化排列后对应的基 矩阵，维数为(M1N1+M2N2)×KL。此时，多频带多视角双雷达ISAR融合成像问题 可转化为一个稀疏表示问题
Figure BDA0002898636170000086
其中，
Figure BDA0002898636170000087
表示雷达1和雷达2观测回波矢量化数据，
Figure BDA0002898636170000088
表示回波数据对应的 基矩阵，a表示成像场景离散网格对应的散射系数幅度矢量。式(7)的稀疏求 解问题可转化为如下的优化问题
Figure BDA0002898636170000091
由于ISAR成像满足稀疏性，故可采用稀疏表示重构方法，利用少量的已 知观测数据
Figure BDA0002898636170000097
求解出散射系数幅度矢量a，再利用全频带全视角对应的基矩阵F 和a即可补全缺失的频带和角度回波数据，得到全频带全视角雷达回波数据， 等效提高发射带宽和观测角度，从而提高成像的距离和方位分辨率。
FLBI算法：
由于(M1N1+M2N2)＜KL，式(8)是一个NP难问题，求解困难。当
Figure BDA0002898636170000098
满足约 束等距性质条件时，可以将带约束的l0范数最小化问题转化为带约束的l1范数最 小化问题：
Figure BDA0002898636170000092
带约束的l1范数最小化问题又可以转化为线性规划问题，然而传统的线性 规划求解方法并不适用于大规模的基矩阵。因此，考虑到噪声存在的情况，将 式(9)的凸优化问题转化为以下正则化形式
Figure BDA0002898636170000093
其中，μ＞0为正则化参数。这样可通过正则项||a||1控制估计值
Figure BDA0002898636170000094
的稀疏性， 通过保真项
Figure BDA0002898636170000095
控制
Figure BDA0002898636170000096
的误差，能同时保证解的稀疏性和准确性，并较好地 抑制噪声。
算法原理：
考虑到多频带多视角ISAR融合成像时矢量化操作后回波矢量和对应的基 矩阵规模较大，不适合选择涉及到大量矩阵求逆的贝叶斯推理类算法求解，故 需要在保证重构精度的同时选取高效快速的重构算法。Bregman迭代类算法在 求解凸优化问题稀疏解时的具有较好的重构性能及快速收敛性，故本申请选择 该类方法进行求解。利用||a||1的Bregman距离代替||a||1，通过Bregman迭代正则 化求解方法，令凸函数J(a)＝||a||1，则有
Figure BDA0002898636170000101
式中，
Figure RE-GDA0002988244120000102
为Bregman距离，基于凸函数J(·)上的u和v两点的Bregman 距离可定义为
Figure BDA0002898636170000104
其中，向量p为次微分
Figure RE-GDA0002988244120000104
中的一个次梯度，<·>表示内积运算。Bregman 迭代实际上可理解为是残量回代的方法，因此本申请以残量迭代模型进行算法 推导。由于Bregman具备残量回代特征，在迭代时会可能导致停滞现象，故引 入一个参数α(0≤α＜1)来控制残量回代权重，对每次的残差施加“惩罚”，则会 加快收敛速度。此时，式(11)可等价为
Figure BDA0002898636170000107
LBI算法是Bregman迭代正则化的线性简化，令
Figure BDA0002898636170000108
对H(a)进 行线性展开，即在ag处泰勒展开，有
Figure BDA0002898636170000109
其中
Figure BDA00028986361700001010
此时，式(11)可进一步表示为
Figure BDA00028986361700001011
令
Figure BDA00028986361700001012
根据平稳点的次微分条件，有
Figure BDA00028986361700001013
则
Figure BDA00028986361700001014
其中，
Figure BDA0002898636170000111
当a＝ag+1时，
Figure BDA0002898636170000112
带入式(15)有
Figure BDA0002898636170000113
对式(16)利用递推公式，有
Figure BDA0002898636170000114
令
Figure BDA0002898636170000115
则有
Figure BDA0002898636170000116
而由式(17)可得
ag+1＝δ(vg+1-μpg+1) (18)
当J(a)＝||a||1时，式(14)可进一步简化为
Figure BDA0002898636170000117
从式(19)可以看出，目标函数的分量是可分离的。对于向量 w＝(w1,w2,…,wM)T，定义向量阈值算子
Figure RE-GDA0002988244120000118
其中，
Figure BDA0002898636170000119
公式(19)的解为
ag+1＝Tμδ(δvg+1)＝δTμ(vg+1) (21)
故凸优化问题的LBI算法的迭代公式为
Figure BDA0002898636170000121
令
Figure BDA0002898636170000122
r0＝0，则有
Figure BDA0002898636170000123
从而有
Figure BDA0002898636170000124
故式(22)可变为以下形式
Figure BDA0002898636170000125
虽然加入权重系数α可在一定程度上加快算法的收敛速度，但若能优化基 矩阵
Figure BDA0002898636170000126
条件数，则可进一步提高收敛速度。基矩阵
Figure BDA0002898636170000127
的条件数定义为：
Figure BDA0002898636170000128
其中，
Figure BDA0002898636170000129
和
Figure BDA00028986361700001210
分别为
Figure BDA00028986361700001211
的最大和最小特征值。从式(24)可 知，条件数不小于1。由于条件数越小，收敛速度越快，且
Figure BDA00028986361700001219
为行满秩矩阵， 故在式(7)两边分别左乘预条件子
Figure BDA00028986361700001212
得到
Figure BDA00028986361700001213
令
Figure BDA00028986361700001214
则有
Figure BDA00028986361700001215
此时，
Figure BDA00028986361700001216
则
Figure BDA00028986361700001217
的条件数为
Figure BDA00028986361700001218
根据式(23)，条件数优化后的FLBI算法的迭代公式为
Figure BDA0002898636170000131
其中，
Figure BDA0002898636170000132
算法实现流程
FLBI算法的基本步骤可总结如下
输入：
Figure BDA0002898636170000133
初始化：a0＝r0＝0，μ＞0，δ＞0，g＝0，
Figure BDA0002898636170000134
第1步：判断是否满足迭代终止条件，若满足则结束迭代转到第2步，若 不满足，则执行：
Figure BDA0002898636170000135
Figure BDA0002898636170000136
Figure BDA0002898636170000137
g＝g+1
第2步：输出迭代结果
Figure BDA0002898636170000138
仿真与实验分析
本申请仿真实验环境为Windows 7 64位操作系统，MATLAB2016A软件平 台，仿真所用计算机主要参数如下：处理器为Intel酷睿i5-6200U，主频为2.30 GHz，内存为4GB。利用ISAR成像二维数据验证算法在多频带多视角ISAR 融合成像中的性能。
仿真1简单模型融合成像性能验证：
为更好地验证本申请算法在实现多频带多视角ISAR融合成像的性能，利 用简单的目标散射点模型进行ISAR融合成像仿真实验。目标中包含12个散射 点，仿真模型如图4a所示，其中，散射点A和B的坐标分别为(0.2,0.2)和(0,0)。 首先产生全频带全视角的ISAR成像回波信号，雷达仿真参数设置如下：雷达 发射线性调频信号，载频为15.26GHz，全频带信号带宽为1.28GHz，采样频率 为1.536GHz，发射脉冲宽度为10us，脉冲重复频率为100Hz，共发射128个脉 冲。在成像时间内，假设目标在距离雷达50Km的高度以5km/s的速度匀速运动，累积转角为7.24°，此时雷达的距离分辨率和方位分辨率分别为0.117m和 0.078m。每个脉冲回波内采样128个距离采样单元的回波数据，经运动补偿和 MTRC校正后作为全频带全视角回波数据，大小为128×128，利用距离-多普勒 (Range Doppler,RD)算法的到的成像结果如图4b所示。
在全频带全视角回波数据中的左上角和右下角分别截取大小为M1×N1和 M2×N2的数据块作为雷达1和雷达2的观测数据，并利用这两段数据进行多频 带多视角融合成像。令M1＝M2＝32，N1＝N2＝32，分别在两部雷达回波数据中加入信 噪比为20dB的高斯白噪声。雷达1和雷达2的RD成像结果分别如图5a和图 5b所示。由于带宽和观测视角有限，单雷达回波的二维分辨率较低，不足以分 辨目标中的散射点。图5c为直接利用两部雷达观测数据进行融合的RD成像结 果，RD融合成像结果的分辨率虽比单雷达结果稍好，但由于频带和视角缺失， 直接FFT进行二维压缩时会引起强烈的副瓣和能量泄露等问题，严重影响了成 像质量。分别利用OMP算法和本申请所提的FLBI算法实现多频带多视角融合 成像，得到的融合成像结果分别如图5d和图5e所示。从图5d可以看出，OMP 算法虽然能较为准确地估计出大多数的目标散射点，但估计的散射点个数和位 置与真实的目标散射点分布相比有偏差，并且容易引入虚假散射点，导致重构 结果不准确。这是因为OMP算法在实现时需要预先设置信号稀疏度，当预设 的稀疏度与目标的散射点个数一致时，重构得到的信号估计效果较为准确，然 而当预设的稀疏度与目标的散射点个数不一致时，重构得到的信号估计误差较大，而在实际过程中，准确地估计散射点个数较为困难，这也制约了算法的应 用。相比之下，本申请所提算法不需要预估目标散射点个数，从图5e可以看出， 利用本申请所提算法能够得到较好的融合成像结果，由于发射信号带宽和观测 视角的等效提高，融合后的二维分辨率也有所改善，能够准确分辨出散射点A 和B。
仿真2复杂模型融合成像性能验证：
为进一步验证算法的融合成像性能，利用波音727飞机的回波数据来验证 算法的有效性。全频带距离向采样点数为64，全视角观测下的ISAR回波脉冲 个数为256。ISAR回波经脉冲压缩和平动补偿，且经过MTRC校正后得到的全视 角全频带距离频域-方位慢时间域二维回波及其RD成像结果分别如图6a和图 6b所示，从其RD成像结果中可以清晰地看出飞机的基本轮廓以及机身部位的 一些细节结构。
在全频带全视角回波数据中的左上角和右下角分别截取大小为M1×N1和 M2×N2的数据块作为雷达1和雷达2的观测数据，并利用这两段数据进行多频 带多视角融合成像。令M1＝M2＝24，N1＝N2＝80，分别在两部雷达回波数据中加入信 噪比为20dB的高斯白噪声。雷达1和雷达2的RD成像结果分别如图7a和图 7b所示，只能从图中大概判断出飞机的机头和机身位置，但其细节信息已经难 以分辨。图7c为直接利用两部雷达观测数据进行融合的RD成像结果，从图中 可以看出机身部位的某些细节结构信息，但由于频带和视角缺失，直接进行FFT 压缩成像时引起了强烈的副瓣和能量泄露等问题，导致飞机轮廓不清晰。分别利用OMP算法和本申请所提的FLBI算法实现多频带多视角融合成像，得到的 融合成像结果分别如图7d和图7e所示。从图7d可以看出，OMP算法虽能重 构出较为干净的飞机基本轮廓，但同时引入了较多的虚假散射点导致机翼和机 身的细节结构不清晰。从图7e可以看出，利用本申请所提算法能够融合得到干 净且清晰的目标轮廓，机头、机翼、机尾以及机身的主体细节结构也清晰可辨， 融合成像效果较好。
为进一步验证算法在不同观测回波数据条件下的融合成像性能，保持观测 数据的SNR为20dB，改变观测数据的带宽和视角大小，分别利用OMP算法 和本申请所提算法实现多视角多频带融合成像。
情形1：令M1＝M2＝16，N1＝N2＝64。雷达1和雷达2的观测二维回波如图8a所 示，利用观测回波数据实现融合成像，利用OMP算法和本申请所提算法得到 的融合成像结果分别如图8b和图8c所示。
情形2：令M1＝M2＝16，N1＝N2＝32。雷达1和雷达2的观测二维回波如图9a所 示，利用观测回波数据实现融合成像，利用OMP算法和本申请所提算法得到 的融合成像结果分别如图9b和图9c所示。
情形3：令M1＝M2＝8，N1＝N2＝64。雷达1和雷达2的观测二维回波如图10a 所示，利用观测回波数据实现融合成像，利用OMP算法和本申请所提算法得 到的融合成像结果分别如图10b和图10c所示。
当M1＝M2＝16，N1＝N2＝64时，从图8b可以看出，利用OMP算法进行融合 成像时，从成像结果中可分辨出飞机的基本形状，但存在一定的虚假散射点破 坏了飞机的整体结构，影响了成像质量。从图8c可以看出，利用FLBI算法进 行融合成像时，从成像结果中可清晰地分辨出目标形状及其细节结构信息，成 像质量较高。随着有效观测回波数据减少，从图9b和图10b可以看出，利用 OMP算法进行融合成像时，成像结果质量严重下降，融合图像中存在大量的虚 假散射点，导致无法从成像结果中分辨出目标的基本形状及其结构信息。相比 而言，从图9c和图10c可以看出，利用FLBI算进行融合成像时，仍能得到较 为清晰的目标图像，可从成像结果中分辨出目标的基本形状。
综上，所述方法在建立多频带多视角融合成像稀疏表示模型的基础上，直 接在复数域利用LBI算法进行求解，同时利用加权残差回代和优化基矩阵条件 数相结合的方式来进一步加快迭代收敛速度，在求解时避免了大量的矩阵求逆 过程，大大提高了成像效率，仿真实验验证了算法的有效性和优越性。

Claims (1)

1.一种多频带多视角ISAR融合成像方法，其特征在于包括如下步骤：
1)对各雷达回波信号进行预处理以及互相干处理，得到相干的距离频域-方位慢时间域回波信号；
2)构造字典矩阵，将回波信号离散化表示；
3)分别将回波数据矢量化拼接，得到待融合的观测信号
Figure FDA0003589276150000011
以及对应的基矩阵
Figure FDA0003589276150000012
4)利用FLBI算法进行迭代求解，得到目标图像矢量估计值
Figure FDA0003589276150000013
5)将一维目标图像矢量估计值
Figure FDA0003589276150000014
转换为二维矩阵
Figure FDA0003589276150000015
Figure FDA0003589276150000018
即为多视角多频带雷达信号融合得到的目标图像；
将所述回波信号离散化表示的方法如下：
设ISAR发射线性调频信号，在远场条件下，经解线频调处理后，雷达回波可表示为：
Figure FDA0003589276150000016
其中，fm为距离频率，tn＝nTr为慢时间函数，n＝0,1,…N-1，N为回波脉冲数，Tr为脉冲重复时间，P为目标散射点个数，ap为散射点p的散射系数，c为电磁传播速率，ΔRp(tn)＝ypcos(Δθ(tn))+xp sin(Δθ(tn))为散射点p到参考点的距离，(xp,yp)为散射点p的坐标，Δθ(tn)为观测累积转角；
由于成像观测时间短，累积转角Δθ(tn)较小，有cos(Δθ(tn))≈1，sin(Δθ(tn))≈Δθn；经运动补偿后，目标运动模型可近似为转台模型，假设匀速转动的角速度为ω，则有Δθ(tn)＝ωtn，则式(1)可近似为：
Figure FDA0003589276150000017
对频率进行离散化采样，令fm＝f0+mΔf，m＝0,1,…M-1，其中，M为频率采样点数，Δf为频率采样间隔；在有限带宽和小角度的情况下，忽略散射点越分辨单元徙动MTRC的影响，雷达回波距离频域可离散表示为：
Figure FDA0003589276150000021
其中，a′p＝apexp(-j4πf0yp/c)；令
Figure FDA0003589276150000022
ISAR成像目标的尺寸较小，有ωm,ωn∈(0,1]，将其离散化，令ωm＝k/K(k＝0,1,…,K-1)，ωn＝l/L(l＝0,1,…,L-1)，且有K≥M，L≥N；此时，式(3)可表示为：
Figure FDA0003589276150000023
其中，akl表示每个像素网格的散射系数幅度；式(4)等价于将目标成像场景离散化，距离向和方位向分别划分为K个和L个网格，一共有K×L个成像像素网格；当某一网格交点的坐标(l,k)上存在等效散射点时，此网格点的幅度akl≠0，反之，当此位置上不存在等效散射点时，则akl＝0；由于目标ISAR图像由有限个散射点组成，仅占成像场景中很小的一部分，即akl中仅有少数幅度为非零，大多数幅度为零，故ISAR图像满足很强的稀疏性；
构造距离向傅里叶变换矩阵FR＝[FR(0) FR(1)…FR(m)…FR(M-1)]T，其中FR(m)＝[exp(-j2πm·0/K) exp(-j2πm·1/K)…exp(-j2πm·(K-1)/K)]，构造方位向傅里叶变换矩阵FA＝[FA(0) FA(1)…FA(n)…FA(N-1)]，其中FA(n)＝[exp(-j2πn·0/L) exp(-j2πn·1/L)…exp(-j2πn·(L-1)/L)]T，则ISAR回波可矩阵表示为
S＝FRΑFA (5)
其中，S为大小为M×N的雷达回波数据，FR为大小为M×K的距离向稀疏基矩阵，FA为大小为L×N的方位向稀疏基矩阵，Α为大小为K×L的散射系数矩阵，该矩阵可代表目标二维ISAR图像，其中第k行l列元素为akl；
得到待融合的观测信号
Figure FDA0003589276150000024
以及对应的基矩阵
Figure FDA0003589276150000025
的方法如下：
两部雷达相近放置，分别工作在不同频带并以不同视角同时对目标进行观测；雷达1发射信号频带为
Figure FDA0003589276150000026
共有M1个频率采样点，对应观测角度为
Figure FDA0003589276150000027
共有N1个角度；雷达2发射信号频带为
Figure FDA0003589276150000028
共有M2个频率采样点，观测角度为
Figure FDA0003589276150000029
共有N2个角度；假设Δf和Δθ分别表示频率采样间隔和角度采样间隔，M和N分别表示全频带全角度的频率采样个数和角度采样个数，则全频带的频率采样数据可表示为fm＝f0+mΔf(m＝0,1,…,M-1)，全角度的角度采样数据可表示为θn＝θ0+nΔθ(n＝0,1,…,N-1)；
将回波二维矩阵按照频率-角度排为一维向量，即将回波S按列矢量化堆叠，可将式(5)的二维信号模型转化为一维矢量形式
Figure FDA0003589276150000031
其中，s＝vect(S)，a＝vect(A)，vect(·)表示将矩阵各列向量依次堆叠为一个矢量的操作，
Figure FDA0003589276150000032
表示矩阵(FA)T和FR的Kronecker乘积，即F为全频带全视角回波矢量对应的基矩阵；
在全频带全视角回波矢量数据s中去除数据缺失的部分即可得到雷达1和雷达2的观测回波矢量数据
Figure FDA00035892761500000310
同时在基矩阵F中也去除缺失数据对应的行，令FA1＝[FA(0) FA(1)…FA(N1-1)]，FA2＝[FA(N-N2) FA(N-N2+1)…FA(N-1)]，FR1＝[FR(0) FR(1)…FR(M1-1)]T，FR2＝[FR(M-M2) FR(M-M2+1)…FR(M-1)]T，则
Figure FDA0003589276150000033
表示雷达1观测回波矢量数据对应的基矩阵，维数为M1N1×KL，
Figure FDA0003589276150000034
表示雷达2观测回波矢量数据对应的基矩阵，维数为M2N2×KL，则
Figure FDA0003589276150000035
表示双雷达观测回波矢量化排列后对应的基矩阵，维数为(M1N1+M2N2)×KL；此时，多频带多视角双雷达ISAR融合成像问题可转化为一个稀疏表示问题
Figure FDA0003589276150000036
其中，
Figure FDA0003589276150000037
表示雷达1和雷达2观测回波矢量化数据，
Figure FDA0003589276150000038
表示回波数据对应的基矩阵，a表示成像场景离散网格对应的散射系数幅度矢量；式(7)的稀疏求解问题可转化为如下的优化问题：
Figure FDA0003589276150000039
由于ISAR成像满足稀疏性，故可采用稀疏表示重构方法，利用少量的已知观测数据
Figure FDA00035892761500000311
求解出散射系数幅度矢量a，再利用全频带全视角对应的基矩阵F和a即可补全缺失的频带和角度回波数据，得到全频带全视角雷达回波数据；
得到目标图像矢量估计值
Figure FDA0003589276150000041
的方法如下：
当
Figure FDA0003589276150000042
满足约束等距性质条件时，可以将带约束的l0范数最小化问题转化为带约束的l1范数最小化问题：
Figure FDA0003589276150000043
带约束的l1范数最小化问题又可以转化为线性规划问题，考虑到噪声存在的情况，将式(9)的凸优化问题转化为以下正则化形式
Figure FDA0003589276150000044
其中，μ＞0为正则化参数；这样可通过正则项||a||1控制估计值
Figure FDA0003589276150000045
的稀疏性，通过保真项
Figure FDA0003589276150000046
控制
Figure FDA0003589276150000047
的误差；
迭代公式的计算方法如下：
利用||a||1的Bregman距离代替||a||1，令凸函数J(a)＝||a||1，则有
Figure FDA0003589276150000048
式中，
Figure FDA0003589276150000049
为Bregman距离，基于凸函数J(·)上的u和v两点的Bregman距离可定义为
Figure FDA00035892761500000410
其中，向量p为次微分
Figure FDA00035892761500000412
中的一个次梯度，<·>表示内积运算；由于Bregman具备残量回代特征，在迭代时会可能导致停滞现象，故引入一个参数α来控制残量回代权重，其中0≤α＜1，对每次的残差施加惩罚，则会加快收敛速度；此时，式(11)可等价为
Figure FDA00035892761500000411
令
Figure FDA0003589276150000051
对H(a)进行线性展开，即在ag处泰勒展开，有
Figure FDA0003589276150000052
其中
Figure FDA0003589276150000053
此时，式(11)可进一步表示为：
Figure FDA0003589276150000054
令
Figure FDA0003589276150000055
根据平稳点的次微分条件，有
Figure FDA0003589276150000056
则
Figure FDA0003589276150000057
其中，
Figure FDA0003589276150000058
当a＝ag+1时，
Figure FDA0003589276150000059
带入式(15)有
Figure FDA00035892761500000510
对式(16)利用递推公式，有
Figure FDA00035892761500000511
令
Figure FDA00035892761500000512
则有
Figure FDA00035892761500000513
而由式(17)可得
ag+1＝δ(vg+1-μpg+1) (18)
当J(a)＝||a||1时，式(14)可进一步简化为
Figure FDA00035892761500000514
从式(19)可以看出，目标函数的分量是可分离的；对于向量w＝(w1,w2,…,wM)T，定义向量阈值算子
Figure FDA0003589276150000061
其中，
Figure FDA0003589276150000062
公式(19)的解为
ag+1＝Tμδ(δvg+1)＝δTμ(vg+1) (21)
故凸优化问题的LBI算法的迭代公式为
Figure FDA0003589276150000063
令
Figure FDA0003589276150000064
r0＝0，则有
Figure FDA0003589276150000065
从而有
Figure FDA0003589276150000066
故式(22)可变为以下形式
Figure FDA0003589276150000067
虽然加入权重系数α可在一定程度上加快算法的收敛速度，但若能优化基矩阵
Figure FDA0003589276150000068
条件数，则可进一步提高收敛速度；
Figure FDA0003589276150000069
其中，
Figure FDA00035892761500000610
和
Figure FDA00035892761500000611
分别为
Figure FDA00035892761500000612
的最大和最小特征值；从式(24)可知，条件数不小于1；由于条件数越小，收敛速度越快，且
Figure FDA00035892761500000613
为行满秩矩阵，故在式(7)两边分别左乘预条件子
Figure FDA00035892761500000614
得到
Figure FDA0003589276150000071
令
Figure FDA0003589276150000072
则有
Figure FDA0003589276150000073
此时，
Figure FDA0003589276150000074
则
Figure FDA0003589276150000075
的条件数为
Figure FDA0003589276150000076
根据式(23)，条件数优化后的FLBI算法的迭代公式为
Figure FDA0003589276150000077
其中，
Figure FDA0003589276150000078
CN202110049776.2A 2021-01-14 2021-01-14 多频带多视角isar融合成像方法 Active CN112859074B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110049776.2A CN112859074B (zh) 2021-01-14 2021-01-14 多频带多视角isar融合成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110049776.2A CN112859074B (zh) 2021-01-14 2021-01-14 多频带多视角isar融合成像方法

Publications (2)

Publication Number Publication Date
CN112859074A CN112859074A (zh) 2021-05-28
CN112859074B true CN112859074B (zh) 2022-07-19

Family

ID=76006152

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110049776.2A Active CN112859074B (zh) 2021-01-14 2021-01-14 多频带多视角isar融合成像方法

Country Status (1)

Country Link
CN (1) CN112859074B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114488152B (zh) * 2022-04-18 2022-07-01 南京信息工程大学 基于后向投影的高效近场大小尺寸目标isar成像方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108333587A (zh) * 2018-02-12 2018-07-27 电子科技大学 基于分裂Bregman的前视扫描雷达超分辨成像方法
CN109633646A (zh) * 2019-01-21 2019-04-16 中国人民解放军陆军工程大学 一种基于加权l1范数约束的双基地isar成像方法
CN110780298A (zh) * 2019-11-01 2020-02-11 西安电子科技大学 基于变分贝叶斯学习的多基isar融合成像方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9713012B2 (en) * 2015-07-21 2017-07-18 RadComm, Inc. Methods, devices and systems for enabling simultaneous operation of different technology based devices over a shared frequency spectrum
CN107340518B (zh) * 2017-07-19 2019-05-24 电子科技大学 一种用于信号缺失下的isar雷达成像方法
CN109633647B (zh) * 2019-01-21 2022-02-08 中国人民解放军陆军工程大学 一种双基地isar稀疏孔径成像方法
CN109782279A (zh) * 2019-01-21 2019-05-21 中国人民解放军陆军工程大学 一种基于压缩感知的双基地isar成像方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108333587A (zh) * 2018-02-12 2018-07-27 电子科技大学 基于分裂Bregman的前视扫描雷达超分辨成像方法
CN109633646A (zh) * 2019-01-21 2019-04-16 中国人民解放军陆军工程大学 一种基于加权l1范数约束的双基地isar成像方法
CN110780298A (zh) * 2019-11-01 2020-02-11 西安电子科技大学 基于变分贝叶斯学习的多基isar融合成像方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
A Bistatic ISAR Sparse Aperture High Resolution Imaging Algorithm based on ROMP Algorithm;Wenhua Hu et al.;《2019 IEEE 3rd Information Technology, Networking, Electronic and Automation Control Conference (ITNEC)》;20190606;全文 *
一种快速复数线性Bregman迭代算法及其在ISAR成像中的应用;李少东等;《中国科学(信息科学)》;20151231;第45卷(第009期);全文 *
任意稀疏结构的复稀疏信号快速重构算法及其逆合成孔径雷达成像;陈文峰等;《光电子激光》;20150415(第04期);全文 *
低信噪比下二维联合快速超分辨B-ISAR成像方法;陈文峰等;《电子学报》;20180415(第04期);全文 *

Also Published As

Publication number Publication date
CN112859074A (zh) 2021-05-28

Similar Documents

Publication Publication Date Title
CN103713288B (zh) 基于迭代最小化稀疏贝叶斯重构线阵sar成像方法
CN109116311B (zh) 基于知识辅助稀疏迭代协方差估计的杂波抑制方法
CN105137424B (zh) 一种杂波背景下实波束扫描雷达角超分辨方法
CN103744076B (zh) 基于非凸优化的mimo雷达动目标检测方法
CN110244303B (zh) 基于sbl-admm的稀疏孔径isar成像方法
CN102749621B (zh) 一种双基地合成孔径雷达频域成像方法
CN103869311A (zh) 实波束扫描雷达超分辨成像方法
CN104950305A (zh) 一种基于稀疏约束的实波束扫描雷达角超分辨成像方法
CN110726992B (zh) 基于结构稀疏和熵联合约束的sa-isar自聚焦法
CN109613532B (zh) 一种机载雷达实时多普勒波束锐化超分辨成像方法
Zhang et al. Matrix completion for downward-looking 3-D SAR imaging with a random sparse linear array
CN108226891B (zh) 一种扫描雷达回波计算方法
CN105699969A (zh) 基于广义高斯约束的最大后验估计角超分辨成像方法
CN104122549B (zh) 基于反卷积的雷达角超分辨成像方法
CN105137425A (zh) 基于卷积反演原理的扫描雷达前视角超分辨方法
CN106291543A (zh) 一种运动平台扫描雷达超分辨成像方法
CN109444882B (zh) 基于变斜视椭圆波束同步模型的双站sar成像方法
CN110879391B (zh) 基于电磁仿真和弹载回波仿真的雷达图像数据集制作方法
Gao et al. Perception through 2d-mimo fmcw automotive radar under adverse weather
CN112859074B (zh) 多频带多视角isar融合成像方法
CN103076608B (zh) 轮廓增强的聚束式合成孔径雷达成像方法
CN116359921A (zh) 基于加速轨迹双基前视合成孔径雷达的快速时域成像方法
CN116184343A (zh) 基于相控阵雷达的三维空间蜂群目标检测及信息估计方法
Thammakhoune et al. Moving target imaging for synthetic aperture radar via RPCA
CN115480245A (zh) 一种基于二维超分辨算法的机载雷达前视三维成像技术

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant