CN114638874B - 基于因式分解和isea的空间目标三维重构方法 - Google Patents

基于因式分解和isea的空间目标三维重构方法 Download PDF

Info

Publication number
CN114638874B
CN114638874B CN202210543687.8A CN202210543687A CN114638874B CN 114638874 B CN114638874 B CN 114638874B CN 202210543687 A CN202210543687 A CN 202210543687A CN 114638874 B CN114638874 B CN 114638874B
Authority
CN
China
Prior art keywords
isar image
target
image sequence
isar
scattering
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
CN202210543687.8A
Other languages
English (en)
Other versions
CN114638874A (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.)
Xidian University
Original Assignee
Xidian University
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 Xidian University filed Critical Xidian University
Priority to CN202210543687.8A priority Critical patent/CN114638874B/zh
Publication of CN114638874A publication Critical patent/CN114638874A/zh
Application granted granted Critical
Publication of CN114638874B publication Critical patent/CN114638874B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/50—Depth or shape recovery
    • G06T7/55—Depth or shape recovery from multiple images
    • 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/9004—SAR image acquisition techniques
    • G01S13/9019—Auto-focussing of the SAR signals
    • 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]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • G06T7/66—Analysis of geometric attributes of image moments or centre of gravity
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G06T2207/10044—Radar image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20081—Training; Learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Electromagnetism (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Data Mining & Analysis (AREA)
  • General Health & Medical Sciences (AREA)
  • Biomedical Technology (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Geometry (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于因式分解和ISEA的空间目标三维重构方法，包括：对空间目标大角度ISAR回波数据进行预处理和成像，得到若干帧ISAR图像序列；对ISAR图像序列中的部分关键特征点进行提取得到航迹矩阵，并利用因式分解方法估计投影矩阵；以散射点在ISAR图像序列上的能量积累值为优化的目标函数，利用粒子群优化算法搜索重构空间目标散射中心的三维分布，得到散射中心的初始集合；基于投影矩阵，将初始集合中的散射中心投影到ISAR图像序列中相应的位置，以对初始集合中的错误点进行剔除，从而得到最终的空间目标三维重构散射中心集合。本发明大幅降低了散射中心提取和关联的难度，并克服了现有ISEA算法无法对运动参数未知的空间目标进行重构的问题。

Description

基于因式分解和ISEA的空间目标三维重构方法
技术领域
本发明属于雷达信号处理技术领域，具体涉及一种基于因式分解和ISEA的空间目标三维重构方法。
背景技术
逆合成孔径雷达（Inverse Synthetic Aperture Radar，ISAR）成像技术是对空间目标进行观测的最有效的途径之一。它通过雷达发射一系列宽带电磁脉冲信号对空间目标进行长时间、大角度的持续观测，并对回波信号进行距离向脉冲压缩和方位向相干积累，从而得到观测目标的二维高分辨率图像。但是，通过上述方法得到的二维ISAR图像只是空间目标三维结构在雷达成像平面的投影，无法真实反映空间目标的三维结构，难以满足后续的空间目标分类和识别的要求，因此，对空间目标三维成像方法的研究逐步成为目前ISAR成像领域的热点。
针对ISAR空间目标三维成像问题，现有技术提供了以下几种方法。第一种是联合方位定标和三维重构的方法，首先，将目标的等效转动角速度建模为关于时间的多项式形式，随后，提出松弛约束的因式分解方法对散射中心三维位置进行重构，最后通过投影向量来估计目标等效旋转运动参数，重新对散射中心方位向进行定标处理。该方法通过迭代应用松弛约束因式分解方法和旋转运动参数估计实现方位定标和散射中心三维位置的重构。但是，由于空间目标在微波频段的电磁散射特性具有各向异性，对于具有复杂结构的空间目标，其三维结构很难等效为固定位置的散射中心，使得散射中心提取和关联的基础不复存在。同时，空间目标不同部件之间存在相互遮挡，也给现有散射中心的提取和关联带来诸多挑战。
第二种是基于ISAR图像序列能量积累（ISEA）的三维几何重建方法。通过构造三轴稳定空间目标三维几何体与ISAR图像序列之间的投影矩阵，利用粒子群优化算法优化求解空间目标的三维结构。然而，当目标本身存在未知的运动时，无法通过对目标进行有效的运动建模来获得投影矩阵，使得ISEA方法失效。
发明内容
为了解决现有技术中存在的上述问题，本发明提供了一种基于因式分解和ISEA的空间目标三维重构方法。本发明要解决的技术问题通过以下技术方案实现：
本发明提供了一种基于因式分解和ISEA的空间目标三维重构方法，包括以下步骤：
步骤1：对空间目标大角度ISAR回波数据进行预处理和成像，得到若干帧ISAR图像序列；
步骤2：对所述ISAR图像序列中的部分关键特征点进行提取得到航迹矩阵，并利用因式分解方法估计投影矩阵；
步骤3：以散射中心在ISAR图像序列上的能量积累值为优化的目标函数，利用粒子群优化算法搜索重构空间目标散射中心的三维分布，得到散射中心的初始集合；
步骤4：基于所述投影矩阵，将所述初始集合中的散射中心投影到所述ISAR图像序列中相应的位置，以对初始集合中的错误点进行剔除，从而得到最终的空间目标三维重构散射中心集合。
在本发明的一个实施例中，步骤1包括：
1.1) 将所述空间目标大角度ISAR回波数据划分为F个子孔径；
1.2) 分别对每个子孔径数据进行包络对齐和自聚焦处理，以实现对所述ISAR回波数据的平动补偿；
1.3) 采用RD算法对经过平动补偿的ISAR回波数据进行成像，得到F帧高分辨二维ISAR图像序列。
在本发明的一个实施例中，步骤2包括：
2.1) 从所述ISAR图像序列中提取P个关键特征点，得到距离-多普勒测量矩阵M；
2.2) 对所述距离-多普勒测量矩阵进行奇异值分解；
2.3) 根据距离投影向量和多普勒投影向量引入约束条件，以对所述距离-多普勒测量矩阵的奇异值分解式进行求解，得到投影矩阵。
在本发明的一个实施例中，步骤3包括：
3.1) 计算ISAR图像序列的总能量
Figure 101971DEST_PATH_IMAGE001
，并初始化ISAR图像序列剩余能量
Figure 81428DEST_PATH_IMAGE002
，初始化重构空间目标散射中心的个数上限为
Figure 48116DEST_PATH_IMAGE003
，初始化ISAR图像序列剩余能量与总能量比值
Figure 381008DEST_PATH_IMAGE004
的下限为
Figure 650840DEST_PATH_IMAGE005
、初始化重构空间目标散射中心集合
Figure 332357DEST_PATH_IMAGE006
；
3.2) 构造粒子群优化算法的适应度函数；
3.3) 利用粒子群优化算法搜索重构空间中目标散射中心的三维位置，并根据每个粒子的适应度函数值更新全局最优位置；
3.4) 重复上述步骤3.3)，直至达到最大迭代次数，输出当前全局最优位置
Figure 537074DEST_PATH_IMAGE007
，并更新空间目标的散射中心位置集合
Figure 922925DEST_PATH_IMAGE008
；
3.5) 更新ISAR图像序列剩余能量
Figure 716437DEST_PATH_IMAGE009
为集合
Figure 647484DEST_PATH_IMAGE010
中的每个散射中心在ISAR图像序列上的能量积累值，并在判断
Figure 579975DEST_PATH_IMAGE011
或集合
Figure 254670DEST_PATH_IMAGE010
中的散射中心个数等于
Figure 902689DEST_PATH_IMAGE012
时，输出当前空间目标的散射中心位置集合作为散射中心的初始集合；否则，返回步骤3.3)，继续对散射中心进行搜索。
在本发明的一个实施例中，在步骤3.2)中，所述粒子群优化算法的适应度函数表示为：
Figure 519484DEST_PATH_IMAGE013
其中，
Figure 885743DEST_PATH_IMAGE014
是候选散射中心的位置，
Figure 98550DEST_PATH_IMAGE015
是第
Figure 931901DEST_PATH_IMAGE016
帧ISAR图像；
Figure 391701DEST_PATH_IMAGE017
和
Figure 58306DEST_PATH_IMAGE018
分别为第
Figure 261754DEST_PATH_IMAGE019
帧成像平面的多普勒投影向量和距离投影向量，
Figure 946682DEST_PATH_IMAGE020
和
Figure 845893DEST_PATH_IMAGE021
分别是ISAR图像
Figure 734214DEST_PATH_IMAGE022
的距离分辨单元和方位分辨单元，
Figure 741354DEST_PATH_IMAGE023
和
Figure 280788DEST_PATH_IMAGE024
分别为
Figure 161019DEST_PATH_IMAGE025
的距离点数和方位点数。
在本发明的一个实施例中，所述步骤3.3)包括：
3.3a) 初始化粒子个数
Figure 723588DEST_PATH_IMAGE026
和最大迭代次数
Figure 271768DEST_PATH_IMAGE027
，令迭代次数
Figure 868972DEST_PATH_IMAGE028
，分别初始化每个粒子的初始位置为
Figure 169372DEST_PATH_IMAGE029
，个体局部最优位置为
Figure 219236DEST_PATH_IMAGE030
，个体局部最优适应度为
Figure 381227DEST_PATH_IMAGE031
，
Figure 582007DEST_PATH_IMAGE032
；
3.3b) 找出适应度最大的个体对应的位置作为全局最优位置
Figure 53308DEST_PATH_IMAGE033
，初始化个体飞行速度为
Figure 590469DEST_PATH_IMAGE034
；
3.3c) 迭代次数加1，更新每个个体的飞行速度为：
Figure 556151DEST_PATH_IMAGE035
以及个体位置为：
Figure 924684DEST_PATH_IMAGE036
其中，
Figure 241920DEST_PATH_IMAGE037
为非负的惯性权值参数，
Figure 79426DEST_PATH_IMAGE038
和
Figure 832488DEST_PATH_IMAGE039
分别是正的加速度常数，
Figure 258790DEST_PATH_IMAGE040
和
Figure 557047DEST_PATH_IMAGE041
分别为服从
Figure 68800DEST_PATH_IMAGE042
之间均匀分布的随机数；
3.3d) 基于所述ISAR图像序列中的每帧图像
Figure 628482DEST_PATH_IMAGE043
，计算每个粒子的适应度：
Figure 456760DEST_PATH_IMAGE045
；
3.3e) 判断每个粒子的适应度与个体局部最优适应度的大小；若
Figure 378449DEST_PATH_IMAGE046
，则更新每个粒子的个体局部最优位置
Figure 174235DEST_PATH_IMAGE047
；
3.3f) 找到当前个体局部最优适应度的最大值，判断其是否大于全局最优位置所对应的适应度，若是，则更新全局最优位置为当前个体局部最优适应度的最大值对应的个体位置；否则，转至步骤3.4)。
在本发明的一个实施例中，在步骤3.5)中，集合
Figure 472362DEST_PATH_IMAGE048
中的每个散射中心在ISAR图像序列上的能量积累值的计算方法为：
3.5a) 计算集合
Figure DEST_PATH_IMAGE049
中的点
Figure 610606DEST_PATH_IMAGE050
在每帧ISAR图像中的投影位置
Figure 703196DEST_PATH_IMAGE051
；
3.5b) 将该投影位置
Figure 737011DEST_PATH_IMAGE052
为中心的邻域
Figure 635566DEST_PATH_IMAGE053
内的能量设置为零，以获得每帧ISAR图像对应的残差图像；
3.5c) 计算当前散射中心对应的所有残差图像的能量总和，并将其作为该散射中心在ISAR图像序列上的能量积累值。
本发明的有益效果：
本发明提供的基于因式分解和ISEA的空间目标三维重构方法仅通过对ISAR图像序列上的关键特征点进行提取，得到距离-多普勒测量矩阵，大大降低了实际中的操作难度；然后利用因式分解估计出投影矩阵，最后采用ISEA方法实现空间目标三维重构，大幅降低了散射中心提取和关联的难度；并且由于本发明方法的投影矩阵是通过图像序列分解得到，其无需事先获知空间目标的在轨运动状态，解决了传统ISEA方法无法重构存在未知在轨姿态空间目标的问题，在实际应用中的可行性和鲁棒性更强。
以下将结合附图及实施例对本发明做进一步详细说明。
附图说明
图1是本发明实施例提供的一种基于因式分解和ISEA的空间目标三维重构方法的流程示意图；
图2是本发明实施例提供的空间目标观测模型示意图；
图3是本发明实施例提供的空间目标点模型；
图4是本发明实施例提供的空间目标点模型的部分高分辨二维ISAR图像；其中，4(a)为第16帧成像结果图；4(b)为第68帧成像结果图；
图5是采用本发明的方法和传统的ISEA方法对图3所示的空间目标点模型的重构结果的对比图；其中，5(a)为采用8个关键点的重构结果；5(b)为采用传统ISEA方法重构的结果；5(c)为采用本发明方法重构的结果。
具体实施方式
下面结合具体实施例对本发明做进一步详细的描述，但本发明的实施方式不限于此。
实施例一
为了解决现有技术中存在的散射中心提取和关联难度高以及传统ISEA方法无法重构存在未知在轨姿态空间目标的问题，本发明提出了一种基于因式分解和ISEA的空间目标三维重构方法，通过对空间目标ISAR图像序列上关键特征点的提取，形成距离-多普勒测量矩阵，然后利用因式分解方法估计空间目标三维结构到二维ISAR图像序列的投影矩阵，然后进一步采用基于ISEA的三维重构方法获得空间目标三维结构。
具体的，请参见图1，图1是本发明实施例提供的一种基于因式分解和ISEA的空间目标三维重构方法的流程示意图，其包括以下步骤：
步骤1：对空间目标大角度ISAR回波数据进行预处理和成像，得到若干帧ISAR图像序列。
1.1) 将空间目标大角度ISAR回波数据划分为F个子孔径。
具体的，首先设定子孔径长度
Figure 625388DEST_PATH_IMAGE054
，将雷达接收机获取的大角度长时间的宽带回波划分为F个子块，每个子块被称为子孔径。
1.2) 分别对每个子孔径数据进行包络对齐和自聚焦处理，以实现对所述ISAR回波数据的平动补偿。
具体的，首先对每个子孔径数据采用相邻相关法进行包络对齐操作，以消除目标相对雷达的平动带来的包络偏移。然后通过基于最小熵准则的自聚焦算法，补偿平动引起的初相误差。
1.3) 采用RD(Rang-Doppler，距离多普勒)成像算法对经过平动补偿的ISAR回波数据进行处理，得到F帧高分辨二维ISAR图像序列。
需要说明的是，在对回波数据进行平动补偿之前，还包括：对接收的雷达回波依次进行高速补偿和距离压缩处理。具体实现方式可参考现有相关技术，本实施例对此不作详细描述。
步骤2：对ISAR图像序列中的部分关键特征点进行提取得到航迹矩阵，并利用因式分解方法估计投影矩阵。
2.1) 从ISAR图像序列中提取P个关键特征点，得到距离-多普勒测量矩阵M。
在本实施例中，可以根据待测目标的形状人为的设置少量关键特征点，对这些关键特征点进行手动提取，得到如下距离-多普勒测量矩阵：
Figure 701928DEST_PATH_IMAGE055
其中，
Figure 401201DEST_PATH_IMAGE056
表示第
Figure 837867DEST_PATH_IMAGE057
个特征点在第
Figure 760824DEST_PATH_IMAGE058
帧ISAR图像上的多普勒和距离测量值，
Figure 195216DEST_PATH_IMAGE059
表示图像序列中的帧索引，
Figure 452891DEST_PATH_IMAGE060
表示特征点的索引，
Figure 443981DEST_PATH_IMAGE061
表示特征点的水平坐标
Figure 411324DEST_PATH_IMAGE062
集合，即特征点相对于转台中心的多普勒测量值；
Figure 78935DEST_PATH_IMAGE063
表示特征点的垂直坐标
Figure 309059DEST_PATH_IMAGE064
集合，即特征点相对于转台中心的距离测量值。
2.2) 对距离-多普勒测量矩阵进行奇异值分解。
具体地，根据
Figure 290791DEST_PATH_IMAGE065
的低秩特性，进行奇异值分解
Figure 172028DEST_PATH_IMAGE066
其中，
Figure 761272DEST_PATH_IMAGE067
是
Figure 934152DEST_PATH_IMAGE068
的对角矩阵，
Figure 453995DEST_PATH_IMAGE069
和
Figure 189739DEST_PATH_IMAGE070
是由对应于奇异值
Figure 215464DEST_PATH_IMAGE071
,
Figure 607131DEST_PATH_IMAGE072
和
Figure 995911DEST_PATH_IMAGE073
的特征向量组成的矩阵，所以，
Figure 336894DEST_PATH_IMAGE074
和
Figure 720471DEST_PATH_IMAGE075
构成
Figure 661751DEST_PATH_IMAGE076
2.3) 根据距离投影向量和多普勒投影向量引入约束条件，以对距离-多普勒测量矩阵的奇异值分解式进行求解，得到投影矩阵。
由于对于任意可逆矩阵
Figure 523396DEST_PATH_IMAGE077
，
Figure 984465DEST_PATH_IMAGE078
和
Figure 615908DEST_PATH_IMAGE079
都可以满足
Figure 716588DEST_PATH_IMAGE080
，即上述分解并不唯一。因此引入以下约束条件：
Figure 194974DEST_PATH_IMAGE081
其中，
Figure 963079DEST_PATH_IMAGE082
和
Figure 485196DEST_PATH_IMAGE083
分别为距离投影向量和多普勒投影向量。
令
Figure 341681DEST_PATH_IMAGE084
则约束条件
Figure 358178DEST_PATH_IMAGE085
和
Figure 43106DEST_PATH_IMAGE086
可写成如下矩阵形式：
Figure 673808DEST_PATH_IMAGE087
其中，
Figure 827709DEST_PATH_IMAGE088
以及
Figure 897165DEST_PATH_IMAGE089
由此，可计算得到
Figure 984070DEST_PATH_IMAGE090
基于
Figure 257444DEST_PATH_IMAGE091
便可构造对称矩阵
Figure 882329DEST_PATH_IMAGE092
，在
Figure 240629DEST_PATH_IMAGE092
满足正定矩阵的条件下，对
Figure 696887DEST_PATH_IMAGE093
进行Cholesky分解得到
Figure 872653DEST_PATH_IMAGE094
，从而得到投影矩阵
Figure 735567DEST_PATH_IMAGE095
和所提取的关键特征点的三维坐标：
Figure 87439DEST_PATH_IMAGE096
可以理解的是，在本实施例中，根据不同的关键特征点以及对应的ISAR图像序列中的不同图像帧，可以得到不同的投影矩阵。显然，所选的关键特征点越多，最终得到的三维重建的质量越好。然而，选择的关键特征点过多，或者选取的图像帧数越多，也就意味着算法越复杂。因此，本实施例优选两帧ISAR图像进行特征点提取。
步骤3：以散射中心(散射点)在ISAR图像序列上的能量积累值为优化的目标函数，利用粒子群优化算法搜索重构空间目标散射中心的三维分布，得到散射中心的初始集合。
本实施例主要利用了空间目标真实三维点被投影到ISAR图像序列时，其总是位于ISAR图像能量较大的目标区域的特点，采用PSO算法(粒子群优化算法)搜索重构空间目标散射点三维分布。
具体的，步骤3包括：
3.1) 计算ISAR图像序列的总能量
Figure 211252DEST_PATH_IMAGE097
，并初始化ISAR图像序列剩余能量
Figure 823499DEST_PATH_IMAGE098
，初始化重构空间目标散射中心的个数上限为
Figure 422977DEST_PATH_IMAGE099
，初始化ISAR图像序列剩余能量与总能量比值
Figure 310030DEST_PATH_IMAGE100
的下限为
Figure 288351DEST_PATH_IMAGE101
、初始化重构空间目标散射点集合
Figure 393DEST_PATH_IMAGE102
。
3.2) 构造粒子群优化算法的适应度函数；其表达式为：
Figure 696953DEST_PATH_IMAGE103
其中，
Figure 590960DEST_PATH_IMAGE104
是候选散射中心的位置，
Figure 610738DEST_PATH_IMAGE105
是第
Figure 564787DEST_PATH_IMAGE106
帧ISAR图像；
Figure 748644DEST_PATH_IMAGE107
和
Figure 980429DEST_PATH_IMAGE108
分别为第
Figure 74287DEST_PATH_IMAGE109
帧成像平面的多普勒投影向量和距离投影向量，
Figure 995976DEST_PATH_IMAGE110
和
Figure 791762DEST_PATH_IMAGE111
分别是ISAR图像
Figure 89889DEST_PATH_IMAGE112
的距离分辨单元和方位分辨单元，
Figure 228133DEST_PATH_IMAGE113
和
Figure 133772DEST_PATH_IMAGE114
分别为
Figure 354538DEST_PATH_IMAGE115
的距离点数和方位点数。
3.3) 利用粒子群优化算法搜索重构空间中目标散射中心的三维位置，并根据每个粒子的适应度函数值更新全局最优位置。
在本实施例中，步骤3.3)主要包括：
3.3a) 初始化粒子个数
Figure 987514DEST_PATH_IMAGE116
和最大迭代次数
Figure 55964DEST_PATH_IMAGE117
，令迭代次数
Figure 319455DEST_PATH_IMAGE118
，分别初始化每个粒子初始位置为
Figure 92763DEST_PATH_IMAGE119
，个体局部最优位置为
Figure 467113DEST_PATH_IMAGE120
，个体局部最优适应度为
Figure 390070DEST_PATH_IMAGE121
，
Figure 886779DEST_PATH_IMAGE122
；
3.3b) 找出适应度最大的个体对应的位置作为全局最优位置
Figure 82137DEST_PATH_IMAGE123
，初始化个体飞行速度为
Figure 73227DEST_PATH_IMAGE124
；
3.3c) 迭代次数h加1，并更新每个个体的飞行速度为：
Figure 52289DEST_PATH_IMAGE126
以及个体位置为：
Figure 719900DEST_PATH_IMAGE127
其中，
Figure 950024DEST_PATH_IMAGE128
为非负的惯性权值参数，
Figure 994072DEST_PATH_IMAGE129
和
Figure 812992DEST_PATH_IMAGE130
分别是正的加速度常数，
Figure 402237DEST_PATH_IMAGE131
和
Figure 637433DEST_PATH_IMAGE132
分别为服从
Figure 157276DEST_PATH_IMAGE042
之间均匀分布的随机数；
3.3d) 基于所述ISAR图像序列中的每帧图像
Figure 643753DEST_PATH_IMAGE133
，计算每个粒子的适应度：
Figure 856428DEST_PATH_IMAGE135
；
3.3e) 判断每个粒子的适应度与个体局部最优适应度的大小；若
Figure 310412DEST_PATH_IMAGE136
，则更新每个粒子的个体局部最优位置
Figure 181416DEST_PATH_IMAGE137
；
3.3f) 找到当前个体局部最优适应度的最大值，判断其是否大于全局最优位置所对应的适应度，若是，则更新全局最优位置为当前个体局部最优适应度的最大值对应的个体位置；否则，转至步骤3.4)。
3.4) 重复上述步骤3.3)，直至达到最大迭代次数，输出当前全局最优位置
Figure 977858DEST_PATH_IMAGE138
，并更新空间目标的散射中心位置集合
Figure 361435DEST_PATH_IMAGE139
。
3.5) 更新ISAR图像序列剩余能量
Figure 37136DEST_PATH_IMAGE140
为集合
Figure 164361DEST_PATH_IMAGE141
中的每个散射中心在ISAR图像序列上的能量积累值，并在判断
Figure 625429DEST_PATH_IMAGE142
或集合
Figure 979575DEST_PATH_IMAGE141
中的散射中心个数等于
Figure 158883DEST_PATH_IMAGE143
时，输出当前空间目标的散射中心位置集合作为初始散射中心集合；否则，返回步骤3.3)，继续对散射中心进行搜索。
在步骤3.5)中，集合
Figure 824220DEST_PATH_IMAGE141
中的每个散射中心在ISAR图像序列上的能量积累值的计算方法如下：
3.5a) 计算集合
Figure 389062DEST_PATH_IMAGE141
中的点
Figure 927491DEST_PATH_IMAGE144
在每帧ISAR图像中的投影位置
Figure 515467DEST_PATH_IMAGE145
；
3.5b) 将该投影位置
Figure 38022DEST_PATH_IMAGE146
为中心的邻域
Figure 660634DEST_PATH_IMAGE147
内的能量设置为零，以获得每帧ISAR图像对应的残差图像；
3.5c) 计算当前散射中心对应的所有残差图像的能量总和，并将其作为该散射中心在ISAR图像序列上的能量积累值。
步骤4：基于投影矩阵，将所述初始集合中的散射中心投影到所述ISAR图像序列中相应的位置，以对初始集合中的错误点进行剔除，从而得到最终的空间目标三维重构散射中心集合。
具体地，基于投影矩阵
Figure 104384DEST_PATH_IMAGE148
，将得到的集合
Figure 507553DEST_PATH_IMAGE149
(也即散射中心的初始集合)投影到ISAR图像序列中相应的位置；计算落入目标区域的投影散射点总数
Figure 514692DEST_PATH_IMAGE150
；如果
Figure DEST_PATH_IMAGE151
，则保留当前点，否则删除该点；最后输出空间目标三维重构散射点集合
Figure 994739DEST_PATH_IMAGE152
。
至此，完成基于因式分解和ISEA的空间目标三维重构。
本发明针对在轨运动状态位置的空间目标三维重构问题，提出了一种联合因式分解和ISAR图像序列能量积累的空间目标三维重构方法，仅需对ISAR图像序列上关键特征点的提取，得到距离-多普勒测量矩阵，进而利用因式分解估计出的投影矩阵，采用ISEA方法实现空间目标三维重构，大幅降低了散射中心提取和关联的难度。此外，由于本发明方法的投影矩阵是通过图像序列分解得到，其无需事先获知空间目标的在轨运动状态。因此本发明方法解决了传统ISEA方法无法重构在轨运动未知的空间目标的问题，在实际应用中的可行性和鲁棒性更强。
实施例二
下面通过点目标仿真成像实验进一步说明本发明提出方法的正确性和有效性。
（1）仿真条件
请参见图2和图3，图2是本发明实施例提供的空间目标观测模型示意图；图3是本发明实施例提供的空间目标点模型。实验的关键参数如表1所示。
Figure 796342DEST_PATH_IMAGE153
（2）仿真实验内容及结果分析
2.1采用高分辨率ISAR成像方法，获得了图2所示的点模型的ISAR图像序列。其结果如图4所示。其中，图4(a)为第16帧成像结果图；图4(b)为第68帧成像结果图；
2.2分别提取每一帧ISAR图像中位于太阳能电池板和卫星SAR天线的顶点的8个特征点，得到距离-多普勒测量矩阵；然后利用因式分解方法获得投影矩阵。8个点的重构结果如图5(a)所示，重构结果过于稀疏，无法反映目标的真实三维结构。基于测量矩阵和投影矩阵，通过本发明所提出的方法重构点模型的三维结构，如图5(c)所示。为了进行比较，使用ISEA方法重构同一点模型的三维结构，如图5(b)所示。
显然，由于目标的未知旋转，使得仅由雷达观测信息构成的投影矩阵失效，进而导致图5(b)所示的传统ISEA方法的重构结果严重失真。此外，传统ISEA方法的重构点数为205个，点云结构过于稀疏，不足以描述点模型的三维结构。但是，如图5(c)所示，本发明所提出方法的重构结果有3000个重构点，能够很好地描述仿真点模型的真实三维结构，与图5(b)中所示的结果相比，本发明方法的重构结果更加稠密和准确，这是因为本发明方法通过因式分解方法构造投影矩阵，目标的运动参数被考虑在内。因此，当空间目标以更复杂的方式运动时，本发明方法具有更强的鲁棒性和更高的精度。
此外，为了对重构结果进行定性分析，本实施例定义了投影图像序列相似性（Similarity of Projected Image Sequence，SPIS）来描述重构结果的准确性，它可以表示为
Figure 171960DEST_PATH_IMAGE155
其中，
Figure 45107DEST_PATH_IMAGE156
表示Hadamard乘积。
Figure 189780DEST_PATH_IMAGE157
和
Figure 427863DEST_PATH_IMAGE158
分别表示第f帧重构结果投影图像和第f帧二值化ISAR图像。
Figure 542975DEST_PATH_IMAGE159
和
Figure 704966DEST_PATH_IMAGE160
分别是
Figure 156676DEST_PATH_IMAGE161
和
Figure 627977DEST_PATH_IMAGE162
中所有像素的能量平均值。m、n和f分别是距离单元、方位单元和图像的索引。
通过计算，ISEA方法和本发明提出方法重构结果的SPIS值分别为0.3613和0.9779，这也验证了本发明提出方法的性能优势。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (6)

1.一种基于因式分解和ISEA的空间目标三维重构方法，其特征在于，包括以下步骤：
步骤1：对空间目标大角度ISAR回波数据进行预处理和成像，得到若干帧ISAR图像序列；
步骤2：对所述ISAR图像序列中的部分关键特征点进行提取得到航迹矩阵，并利用因式分解方法估计投影矩阵；具体包括：
2.1)从所述ISAR图像序列中提取P个关键特征点，得到距离-多普勒测量矩阵M；
2.2)对所述距离-多普勒测量矩阵进行奇异值分解；
2.3)根据距离投影向量和多普勒投影向量引入约束条件，以对所述距离-多普勒测量矩阵的奇异值分解式进行求解，得到投影矩阵；
步骤3：以散射中心在ISAR图像序列上的能量积累值为优化的目标函数，利用粒子群优化算法搜索重构空间目标散射中心的三维分布，得到散射中心的初始集合；
步骤4：基于所述投影矩阵，将所述初始集合中的散射中心投影到所述ISAR图像序列中相应的位置，以对初始集合中的错误点进行剔除，从而得到最终的空间目标三维重构散射中心集合。
2.根据权利要求1所述的基于因式分解和ISEA的空间目标三维重构方法，其特征在于，步骤1包括：
1.1)将所述空间目标大角度ISAR回波数据划分为F个子孔径；
1.2)分别对每个子孔径数据进行包络对齐和自聚焦处理，以实现对所述ISAR回波数据的平动补偿；
1.3)采用RD算法对经过平动补偿的ISAR回波数据进行成像，得到F帧高分辨二维ISAR图像序列。
3.根据权利要求1所述的基于因式分解和ISEA的空间目标三维重构方法，其特征在于，步骤3包括：
3.1)计算ISAR图像序列的总能量Etotal，并初始化ISAR图像序列剩余能量Eremain＝Etotal，初始化重构空间目标散射中心的个数上限为Nmax，初始化ISAR图像序列剩余能量与总能量比值
Figure FDA0003738047580000021
的下限为δ、初始化重构空间目标散射中心集合
Figure FDA0003738047580000022
3.2)构造粒子群优化算法的适应度函数；
3.3)利用粒子群优化算法搜索重构空间中目标散射中心的三维位置，并根据每个粒子的适应度函数值更新全局最优位置；
3.4)重复上述步骤3.3)，直至达到最大迭代次数，输出当前全局最优位置
Figure FDA0003738047580000023
并更新空间目标的散射中心位置集合Θ＝Θ∪{popt}；
3.5)更新ISAR图像序列剩余能量Eremain为集合Θ中的每个散射中心在ISAR图像序列上的能量积累值，并在判断
Figure FDA0003738047580000024
或集合Θ中的散射中心个数等于Nmax时，输出当前空间目标的散射中心位置集合作为散射中心的初始集合；否则，返回步骤3.3)，继续对散射中心进行搜索。
4.根据权利要求3所述的基于因式分解和ISEA的空间目标三维重构方法，其特征在于，在步骤3.2)中，所述粒子群优化算法的适应度函数表示为：
Figure FDA0003738047580000025
其中，p＝[x,y,z]T是候选散射中心的位置，If是第f帧ISAR图像；if和jf分别为第f帧成像平面的多普勒投影向量和距离投影向量，Δr和Δf分别是ISAR图像If的距离分辨单元和方位分辨单元，Mr和Ma分别为If的距离点数和方位点数。
5.根据权利要求3所述的基于因式分解和ISEA的空间目标三维重构方法，其特征在于，所述步骤3.3)包括：
3.3a)初始化粒子个数Υ和最大迭代次数H，令迭代次数h＝1，分别初始化每个粒子的初始位置为
Figure FDA0003738047580000031
个体局部最优位置为
Figure FDA0003738047580000032
个体局部最优适应度为
Figure FDA0003738047580000033
3.3b)找出适应度最大的个体对应的位置作为全局最优位置pgbest＝(xgbest,ygbest,zgbest)T，初始化个体飞行速度为
Figure FDA0003738047580000034
3.3c)迭代次数加1，更新每个个体的飞行速度为：
Figure FDA0003738047580000035
以及个体位置为：
Figure FDA0003738047580000036
其中，α为非负的惯性权值参数，c1和c2分别是正的加速度常数，r1和r2分别为服从[0,1]之间均匀分布的随机数；
3.3d)基于所述ISAR图像序列中的每帧图像If，计算每个粒子的适应度：
Figure FDA0003738047580000037
3.3e)判断每个粒子的适应度与个体局部最优适应度的大小；若
Figure FDA0003738047580000038
则更新每个粒子的个体局部最优位置
Figure FDA0003738047580000039
3.3f)找到当前个体局部最优适应度的最大值，判断其是否大于全局最优位置所对应的适应度，若是，则更新全局最优位置为当前个体局部最优适应度的最大值对应的个体位置；否则，转至步骤3.4)。
6.根据权利要求3所述的基于因式分解和ISEA的空间目标三维重构方法，其特征在于，在步骤3.5)中，集合Θ中的每个散射中心在ISAR图像序列上的能量积累值的计算方法为：
3.5a)计算集合Θ中的点popt在每帧ISAR图像中的投影位置
Figure FDA0003738047580000041
3.5b)将该投影位置
Figure FDA0003738047580000042
为中心的邻域ε内的能量设置为零，以获得每帧ISAR图像对应的残差图像；
3.5c)计算当前散射中心对应的所有残差图像的能量总和，并将其作为该散射中心在ISAR图像序列上的能量积累值。
CN202210543687.8A 2022-05-19 2022-05-19 基于因式分解和isea的空间目标三维重构方法 Active CN114638874B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210543687.8A CN114638874B (zh) 2022-05-19 2022-05-19 基于因式分解和isea的空间目标三维重构方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210543687.8A CN114638874B (zh) 2022-05-19 2022-05-19 基于因式分解和isea的空间目标三维重构方法

Publications (2)

Publication Number Publication Date
CN114638874A CN114638874A (zh) 2022-06-17
CN114638874B true CN114638874B (zh) 2022-09-16

Family

ID=81953047

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210543687.8A Active CN114638874B (zh) 2022-05-19 2022-05-19 基于因式分解和isea的空间目标三维重构方法

Country Status (1)

Country Link
CN (1) CN114638874B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115169229B (zh) * 2022-06-29 2023-04-18 广州海洋地质调查局 一种海底底质类型的划分方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9702971B2 (en) * 2014-03-17 2017-07-11 Raytheon Company High-availability ISAR image formation
CN111157985B (zh) * 2019-11-15 2023-04-21 西安电子科技大学 基于多站一维距离像序列的空间刚体目标三维重构方法
CN111208513B (zh) * 2020-01-15 2023-03-31 西安电子科技大学 空间目标isar图像序列能量反向投影与三维重构方法
CN111830504B (zh) * 2020-07-23 2023-11-24 中山大学 基于序贯融合因式分解的序列isar三维成像方法
CN112782695B (zh) * 2021-01-27 2023-05-30 西安电子科技大学 基于isar图像和参数优化的卫星姿态和尺寸估计方法
CN113640791B (zh) * 2021-06-09 2023-12-26 西安电子科技大学 一种基于距离和瞬时速度的空间目标三维姿态重构方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
A New 3-D Geometry Reconstruction Method of Space Target Utilizing the Scatterer Energy Accumulation of ISAR Image Sequence;Lei Liu et al.;《 IEEE Transactions on Geoscience and Remote Sensing》;20200429;8345-8357 *
逆合成孔径雷达二维及三维成像方法研究;刘磊;《中国优秀博士学位论文全文数据库 (信息科技辑)》;20170215;96-111 *

Also Published As

Publication number Publication date
CN114638874A (zh) 2022-06-17

Similar Documents

Publication Publication Date Title
CN111208513B (zh) 空间目标isar图像序列能量反向投影与三维重构方法
CN110058237B (zh) 面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法
CN111948654B (zh) 机载层析sar三维点云生成方法
CN102393518B (zh) 一种适用于大斜视角的机载sar成像方法
CN110148165B (zh) 一种基于粒子群优化的三维干涉isar图像配准方法
CN111157985B (zh) 基于多站一维距离像序列的空间刚体目标三维重构方法
CN106501865B (zh) 一种边缘嵌套加权的稀疏成像方法
CN105699969A (zh) 基于广义高斯约束的最大后验估计角超分辨成像方法
CN112099011A (zh) 基于focuss算法的全息sar子孔径三维重建方法
CN112415515B (zh) 一种机载圆迹sar对不同高度目标分离的方法
CN110133682A (zh) 星载全方位sar自适应目标三维重建方法
CN110018474A (zh) 基于地球同步轨道合成孔径雷达层析技术的三维成像方法
Han et al. Efficient 3D image reconstruction of airborne TomoSAR based on back projection and improved adaptive ISTA
CN114638874B (zh) 基于因式分解和isea的空间目标三维重构方法
CN114660606B (zh) 低信噪比isar图像序列匹配搜索的空间目标姿态反演方法
CN110596706B (zh) 一种基于三维图像域投射变换的雷达散射截面积外推方法
CN113608218B (zh) 一种基于后向投影原理的频域干涉相位稀疏重构方法
CN107678009B (zh) 干涉处理艇载雷达阵列形变误差补偿和目标探测方法
CN109959933A (zh) 一种基于压缩感知的多基线圆迹合成孔径雷达成像方法
CN114067064A (zh) 基于多视角雷达图像的目标三维重建方法
CN112684446B (zh) 基于最小熵准则的Bi-ISAR横向定标与畸变校正方法
Liu et al. Analysis of Deep Learning 3-D Imaging Methods Based on UAV SAR
Feng et al. Multiview isar imaging for complex targets based on improved sbr scattering model
CN112946601B (zh) 基于Gauss-Seidel的高效分布式目标相位优化方法
CN117932195B (zh) 一种星载sar图像在轨定位迭代初值计算方法

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