CN112859074A - 多频带多视角isar融合成像方法 - Google Patents
多频带多视角isar融合成像方法 Download PDFInfo
- Publication number
- CN112859074A CN112859074A CN202110049776.2A CN202110049776A CN112859074A CN 112859074 A CN112859074 A CN 112859074A CN 202110049776 A CN202110049776 A CN 202110049776A CN 112859074 A CN112859074 A CN 112859074A
- Authority
- CN
- China
- Prior art keywords
- radar
- echo
- isar
- matrix
- imaging
- 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.)
- Granted
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 121
- 230000004927 fusion Effects 0.000 title claims abstract description 80
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 76
- 239000011159 matrix material Substances 0.000 claims abstract description 55
- 238000000034 method Methods 0.000 claims abstract description 41
- 238000005457 optimization Methods 0.000 claims abstract description 12
- 238000006467 substitution reaction Methods 0.000 claims abstract description 9
- 238000005070 sampling Methods 0.000 claims description 18
- 238000002592 echocardiography Methods 0.000 claims description 9
- 230000000007 visual effect Effects 0.000 claims description 9
- 230000008569 process Effects 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 8
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 0.000 claims description 4
- 238000006243 chemical reaction Methods 0.000 claims description 4
- 230000001427 coherent effect Effects 0.000 claims description 4
- 238000007781 pre-processing Methods 0.000 claims description 4
- 230000005012 migration Effects 0.000 claims description 3
- 238000013508 migration Methods 0.000 claims description 3
- 230000002452 interceptive effect Effects 0.000 claims description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 2
- 238000004088 simulation Methods 0.000 abstract description 12
- 238000003672 processing method Methods 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 5
- 230000006835 compression Effects 0.000 description 3
- 238000007906 compression Methods 0.000 description 3
- 230000006870 function Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000012937 correction Methods 0.000 description 2
- 230000008878 coupling Effects 0.000 description 2
- 238000010168 coupling process Methods 0.000 description 2
- 238000005859 coupling reaction Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 241000282465 Canis Species 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 238000007499 fusion processing Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000013139 quantization Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 230000001502 supplementing effect Effects 0.000 description 1
- 238000011426 transformation method Methods 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
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]
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种多频带多视角ISAR融合成像方法,涉及图像处理方法技术领域。所述方法首先,建立多频带多视角ISAR融合成像的矢量化稀疏表示模型,为了减少求解时的运算复杂度,提高成像效率,在Bregman迭代算法的基础上,结合加权残差回代和基矩阵条件数优化进一步提高算法的迭代收敛速度。仿真实验验证了算法在实现一维稀疏信号重构以及多频带多视角ISAR融合成像时有效性和优越性。
Description
技术领域
本发明涉及图像处理方法技术领域,尤其涉及一种多频带多视角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)构造字典矩阵,将回波信号离散化表示;
采用上述技术方案所产生的有益效果在于:所述方法在建立多频带多视角 融合成像稀疏表示模型的基础上,直接在复数域利用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)构造字典矩阵,将回波信号离散化表示;
下面结合具体步骤对以上方法进行详细说明:
多频带多视角ISAR融合成像稀疏表示模型:
在单部雷达***观测条件下,由于有限的发射信号带宽和观测角度,导致 其距离分辨率和方位分辨率受到限制。而在多雷达***观测条件下,若能利用 不同工作频带和不同观测视角的多个雷达同时对目标进行观测,再利用信号处 理技术将观测数据融合成一个更大带宽和更大视角的回波信号,则可同时提高 雷达成像的二维分辨率。值得注意的是,为实现融合成像,多部雷达观测目标 同一散射中心的信息不能差距太大,因此雷达间的观测视角差距不能太大。在 高频段,频带较窄、观测角度较小的情况下,散射点的复幅度可认为是常数, 即不同雷达观测到同一目标的散射响应特性可近似为理想散射点模型。假设所 讨论的多雷达信号已经过预处理和非相干量补偿,实现了信号的完全相关,本 申请以两个不同频带不同视角的雷达信号为例进行ISAR融合成像分析。
单雷达***ISAR成像模型:
假设ISAR发射线性调频信号,在远场条件下,经解线频调处理后,雷达 回波可表示为
其中,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)可近似为
对频率进行离散化采样,令fm=f0+mΔf,m=0,1,…M-1,其中,M为频率采样 点数,Δf为频率采样间隔。在有限带宽和小角度的情况下,忽略散射点越分辨 单元徙动(Migration Through Resolution Cells,MTRC)的影响,雷达回波距离频 域可离散表示为
其中,a′p=apexp(-j4πf0yp/c)。若存在MTRC,可利先利用Keystone变换法重 新定义一个慢时间,以校正由快时间频域-慢时间耦合项引起的越距离单元徙动 校正,再构造一个相位补偿项以校正越多普勒单元徙动。
此时,式(3)可表示为
其中,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发 射信号频带为共有M1个频率采样点,对应观测角度为共有N1个角度;雷达2发射信号频带为共有M2个频率采样点, 观测角度为共有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)的二维信 号模型转化为一维矢量形式
其中,s=vect(S),a=vect(A),vect(·)表示将矩阵各列向量依次堆叠为一个矢量的操作,表示矩阵(FA)T和FR的Kronecker乘积,即F为全频带全 视角回波矢量对应的基矩阵。基于矢量化处理的多频带多视角双雷达观测数据 融合示意图如图2所示。其中红色矩形部分表示与雷达1的N1个观测角度对应 的基,蓝色矩形部分表示与雷达2的N2个观测角度对应的基。
在全频带全视角回波矢量数据s中去除数据缺失的部分即可得到雷达1和 雷达2的观测回波矢量数据同时在基矩阵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,则 表示雷达1观测回波矢量数据对应的基矩阵,维数为M1N1×KL,对 应图2中去除数据缺失部分后的上侧矩形块部分,表示雷达2观 测回波矢量数据对应的基矩阵,维数为M2N2×KL,对应图2中去除数据缺失部 分后的下侧矩形块部分,则表示双雷达观测回波矢量化排列后对应的基 矩阵,维数为(M1N1+M2N2)×KL。此时,多频带多视角双雷达ISAR融合成像问题 可转化为一个稀疏表示问题
由于ISAR成像满足稀疏性,故可采用稀疏表示重构方法,利用少量的已 知观测数据求解出散射系数幅度矢量a,再利用全频带全视角对应的基矩阵F 和a即可补全缺失的频带和角度回波数据,得到全频带全视角雷达回波数据, 等效提高发射带宽和观测角度,从而提高成像的距离和方位分辨率。
FLBI算法:
带约束的l1范数最小化问题又可以转化为线性规划问题,然而传统的线性 规划求解方法并不适用于大规模的基矩阵。因此,考虑到噪声存在的情况,将 式(9)的凸优化问题转化为以下正则化形式
算法原理:
考虑到多频带多视角ISAR融合成像时矢量化操作后回波矢量和对应的基 矩阵规模较大,不适合选择涉及到大量矩阵求逆的贝叶斯推理类算法求解,故 需要在保证重构精度的同时选取高效快速的重构算法。Bregman迭代类算法在 求解凸优化问题稀疏解时的具有较好的重构性能及快速收敛性,故本申请选择 该类方法进行求解。利用||a||1的Bregman距离代替||a||1,通过Bregman迭代正则 化求解方法,令凸函数J(a)=||a||1,则有
其中,向量p为次微分中的一个次梯度,<·>表示内积运算。Bregman 迭代实际上可理解为是残量回代的方法,因此本申请以残量迭代模型进行算法 推导。由于Bregman具备残量回代特征,在迭代时会可能导致停滞现象,故引 入一个参数α(0≤α<1)来控制残量回代权重,对每次的残差施加“惩罚”,则会 加快收敛速度。此时,式(11)可等价为
对式(16)利用递推公式,有
ag+1=δ(vg+1-μpg+1) (18)
当J(a)=||a||1时,式(14)可进一步简化为
从式(19)可以看出,目标函数的分量是可分离的。对于向量 w=(w1,w2,…,wM)T,定义向量阈值算子
ag+1=Tμδ(δvg+1)=δTμ(vg+1) (21)
故凸优化问题的LBI算法的迭代公式为
根据式(23),条件数优化后的FLBI算法的迭代公式为
算法实现流程
FLBI算法的基本步骤可总结如下
第1步:判断是否满足迭代终止条件,若满足则结束迭代转到第2步,若 不满足,则执行:
g=g+1
仿真与实验分析
本申请仿真实验环境为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 (5)
2.如权利要求1所述的多频带多视角ISAR融合成像方法,其特征在于,将所述回波信号离散化表示的方法如下:
设ISAR发射线性调频信号,在远场条件下,经解线频调处理后,雷达回波可表示为:
其中,fm为距离频率,tn=nTr为慢时间函数,n=0,1,…N-1,N为回波脉冲数,Tr为脉冲重复时间,P为目标散射点个数,ap为散射点p的散射系数,c为电磁传播速率,ΔRp(tn)=ypcos(Δθ(tn))+xpsin(Δθ(tn))为散射点p到参考点的距离,(xp,yp)为散射点p的坐标,Δθ(tn)为观测累积转角;
由于成像观测时间短,累积转角Δθ(tn)较小,有cos(Δθ(tn))≈1,sin(Δθ(tn))≈Δθn。经运动补偿后,目标运动模型可近似为转台模型,假设匀速转动的角速度为ω,则有Δθ(tn)=ωtn,则式(1)可近似为:
对频率进行离散化采样,令fm=f0+mΔf,m=0,1,…M-1,其中,M为频率采样点数,Δf为频率采样间隔;在有限带宽和小角度的情况下,忽略散射点越分辨单元徙动MTRC的影响,雷达回波距离频域可离散表示为:
其中,a′p=apexp(-j4πf0yp/c);令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)可表示为:
其中,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。
两部雷达相近放置,分别工作在不同频带并以不同视角同时对目标进行观测;雷达1发射信号频带为共有M1个频率采样点,对应观测角度为共有N1个角度;雷达2发射信号频带为共有M2个频率采样点,观测角度为共有N2个角度;假设Δf和Δθ分别表示频率采样间隔和角度采样间隔,M和N分别表示全频带全角度的频率采样个数和角度采样个数,则全频带的频率采样数据可表示为fm=f0+mΔf(m=0,1,…,M-1),全角度的角度采样数据可表示为θn=θ0+nΔθ(n=0,1,…,N-1);
将回波二维矩阵按照频率-角度排为一维向量,即将回波S按列矢量化堆叠,可将式(5)的二维信号模型转化为一维矢量形式
在全频带全视角回波矢量数据s中去除数据缺失的部分即可得到雷达1和雷达2的观测回波矢量数据同时在基矩阵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,则表示雷达1观测回波矢量数据对应的基矩阵,维数为M1N1×KL,表示雷达2观测回波矢量数据对应的基矩阵,维数为M2N2×KL,则表示双雷达观测回波矢量化排列后对应的基矩阵,维数为(M1N1+M2N2)×KL;此时,多频带多视角双雷达ISAR融合成像问题可转化为一个稀疏表示问题
5.如权利要求4所述的多频带多视角ISAR融合成像方法,其特征在于,迭代公式的计算方法如下:
利用||a||1的Bregman距离代替||a||1,令凸函数J(a)=||a||1,则有
其中,向量p为次微分中的一个次梯度,<·>表示内积运算;由于Bregman具备残量回代特征,在迭代时会可能导致停滞现象,故引入一个参数α(0≤α<1)来控制残量回代权重,对每次的残差施加惩罚,则会加快收敛速度;此时,式(11)可等价为
对式(16)利用递推公式,有
ag+1=δ(vg+1-μpg+1) (18)
当J(a)=||a||1时,式(14)可进一步简化为
从式(19)可以看出,目标函数的分量是可分离的;对于向量w=(w1,w2,…,wM)T,定义向量阈值算子
ag+1=Tμδ(δvg+1)=δTμ(vg+1) (21)
故凸优化问题的LBI算法的迭代公式为
根据式(23),条件数优化后的FLBI算法的迭代公式为
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 true CN112859074A (zh) | 2021-05-28 |
CN112859074B 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) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114488152A (zh) * | 2022-04-18 | 2022-05-13 | 南京信息工程大学 | 基于后向投影的高效近场大小尺寸目标isar成像方法 |
CN115453527A (zh) * | 2022-08-02 | 2022-12-09 | 西安电子科技大学 | 一种周期性分段观测isar高分辨成像方法 |
CN117538940A (zh) * | 2023-10-27 | 2024-02-09 | 西安电子科技大学 | 基于vna的步进频率探地雷达***及工作方法 |
CN118191756A (zh) * | 2024-05-20 | 2024-06-14 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种雷达信号检测方法及*** |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20170026847A1 (en) * | 2015-07-21 | 2017-01-26 | RadComm, Inc. | Methods, devices and systems for enabling simultaneous operation of different technology based devices over a shared frequency spectrum |
CN107340518A (zh) * | 2017-07-19 | 2017-11-10 | 电子科技大学 | 一种用于信号缺失下的isar雷达成像方法 |
CN108333587A (zh) * | 2018-02-12 | 2018-07-27 | 电子科技大学 | 基于***Bregman的前视扫描雷达超分辨成像方法 |
CN109633647A (zh) * | 2019-01-21 | 2019-04-16 | 中国人民解放军陆军工程大学 | 一种双基地isar稀疏孔径成像方法 |
CN109633646A (zh) * | 2019-01-21 | 2019-04-16 | 中国人民解放军陆军工程大学 | 一种基于加权l1范数约束的双基地isar成像方法 |
CN109782279A (zh) * | 2019-01-21 | 2019-05-21 | 中国人民解放军陆军工程大学 | 一种基于压缩感知的双基地isar成像方法 |
CN110780298A (zh) * | 2019-11-01 | 2020-02-11 | 西安电子科技大学 | 基于变分贝叶斯学习的多基isar融合成像方法 |
-
2021
- 2021-01-14 CN CN202110049776.2A patent/CN112859074B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20170026847A1 (en) * | 2015-07-21 | 2017-01-26 | RadComm, Inc. | Methods, devices and systems for enabling simultaneous operation of different technology based devices over a shared frequency spectrum |
CN107340518A (zh) * | 2017-07-19 | 2017-11-10 | 电子科技大学 | 一种用于信号缺失下的isar雷达成像方法 |
CN108333587A (zh) * | 2018-02-12 | 2018-07-27 | 电子科技大学 | 基于***Bregman的前视扫描雷达超分辨成像方法 |
CN109633647A (zh) * | 2019-01-21 | 2019-04-16 | 中国人民解放军陆军工程大学 | 一种双基地isar稀疏孔径成像方法 |
CN109633646A (zh) * | 2019-01-21 | 2019-04-16 | 中国人民解放军陆军工程大学 | 一种基于加权l1范数约束的双基地isar成像方法 |
CN109782279A (zh) * | 2019-01-21 | 2019-05-21 | 中国人民解放军陆军工程大学 | 一种基于压缩感知的双基地isar成像方法 |
CN110780298A (zh) * | 2019-11-01 | 2020-02-11 | 西安电子科技大学 | 基于变分贝叶斯学习的多基isar融合成像方法 |
Non-Patent Citations (4)
Title |
---|
WENHUA HU ET AL.: "A Bistatic ISAR Sparse Aperture High Resolution Imaging Algorithm based on ROMP Algorithm", 《2019 IEEE 3RD INFORMATION TECHNOLOGY, NETWORKING, ELECTRONIC AND AUTOMATION CONTROL CONFERENCE (ITNEC)》 * |
李少东等: "一种快速复数线性Bregman迭代算法及其在ISAR成像中的应用", 《中国科学(信息科学)》 * |
陈文峰等: "任意稀疏结构的复稀疏信号快速重构算法及其逆合成孔径雷达成像", 《光电子激光》 * |
陈文峰等: "低信噪比下二维联合快速超分辨B-ISAR成像方法", 《电子学报》 * |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114488152A (zh) * | 2022-04-18 | 2022-05-13 | 南京信息工程大学 | 基于后向投影的高效近场大小尺寸目标isar成像方法 |
CN115453527A (zh) * | 2022-08-02 | 2022-12-09 | 西安电子科技大学 | 一种周期性分段观测isar高分辨成像方法 |
CN117538940A (zh) * | 2023-10-27 | 2024-02-09 | 西安电子科技大学 | 基于vna的步进频率探地雷达***及工作方法 |
CN118191756A (zh) * | 2024-05-20 | 2024-06-14 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种雷达信号检测方法及*** |
Also Published As
Publication number | Publication date |
---|---|
CN112859074B (zh) | 2022-07-19 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112859074B (zh) | 多频带多视角isar融合成像方法 | |
CN103713288B (zh) | 基于迭代最小化稀疏贝叶斯重构线阵sar成像方法 | |
CN103487802B (zh) | 扫描雷达角超分辨成像方法 | |
CN103869311B (zh) | 实波束扫描雷达超分辨成像方法 | |
CN109116311B (zh) | 基于知识辅助稀疏迭代协方差估计的杂波抑制方法 | |
Rao et al. | Adaptive sparse recovery by parametric weighted L $ _ {1} $ minimization for ISAR imaging of uniformly rotating targets | |
CN104950306B (zh) | 一种海杂波背景下前视海面目标角超分辨成像方法 | |
CN105652273B (zh) | 一种基于混合匹配追踪算法的mimo雷达稀疏成像算法 | |
CN103487803B (zh) | 迭代压缩模式下机载扫描雷达成像方法 | |
CN102323583B (zh) | 一种超分辨线阵三维合成孔径雷达成像方法 | |
CN103744076B (zh) | 基于非凸优化的mimo雷达动目标检测方法 | |
CN105137424B (zh) | 一种杂波背景下实波束扫描雷达角超分辨方法 | |
CN104950305A (zh) | 一种基于稀疏约束的实波束扫描雷达角超分辨成像方法 | |
CN110726992B (zh) | 基于结构稀疏和熵联合约束的sa-isar自聚焦法 | |
Zhang et al. | Matrix completion for downward-looking 3-D SAR imaging with a random sparse linear array | |
CN104122549B (zh) | 基于反卷积的雷达角超分辨成像方法 | |
CN108226891B (zh) | 一种扫描雷达回波计算方法 | |
CN105652271B (zh) | 一种增广拉格朗日实波束雷达角超分辨处理方法 | |
CN105137425A (zh) | 基于卷积反演原理的扫描雷达前视角超分辨方法 | |
CN106918810B (zh) | 一种存在阵元幅相误差时的微波关联成像方法 | |
CN109597075B (zh) | 一种基于稀疏阵列的成像方法及成像装置 | |
CN109444882B (zh) | 基于变斜视椭圆波束同步模型的双站sar成像方法 | |
CN105699969A (zh) | 基于广义高斯约束的最大后验估计角超分辨成像方法 | |
CN106291543A (zh) | 一种运动平台扫描雷达超分辨成像方法 | |
CN109613532A (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 |