CN115453528A - 基于快速sbl算法实现分段观测isar高分辨成像方法及装置 - Google Patents

基于快速sbl算法实现分段观测isar高分辨成像方法及装置 Download PDF

Info

Publication number
CN115453528A
CN115453528A CN202210944216.8A CN202210944216A CN115453528A CN 115453528 A CN115453528 A CN 115453528A CN 202210944216 A CN202210944216 A CN 202210944216A CN 115453528 A CN115453528 A CN 115453528A
Authority
CN
China
Prior art keywords
matrix
value
algorithm
inverse
segmented
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.)
Pending
Application number
CN202210944216.8A
Other languages
English (en)
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 CN202210944216.8A priority Critical patent/CN115453528A/zh
Publication of CN115453528A publication Critical patent/CN115453528A/zh
Pending legal-status Critical Current

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/9004—SAR image acquisition techniques
    • 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)
  • Image Analysis (AREA)

Abstract

本发明涉及一种基于快速SBL算法实现分段观测ISAR高分辨成像方法，包括：对分段观测数据进行建模，得到分段观测数据模型；基于分段观测数据模型下的傅里叶字典矩阵，结合信号精度值和噪声精度值，采用快速傅里叶变换计算稀疏贝叶斯学习算法进行迭代过程中的待求逆矩阵；基于待求逆矩阵，利用Levinson‑Durbin迭代算法求解逆矩阵的G‑S分解因子；结合逆矩阵的G‑S分解因子，利用快速傅里叶算法求解稀疏贝叶斯学习算法迭代过程中后验分布的均值和协方差的对角元素；当判断均值不满足收敛条件时，计算新的信号精度值和新的噪声精度值，并返回步骤S2；当判断均值满足收敛条件时，输出均值。该成像方法在保证成像结果准确性的基础上减小计算复杂度，提高了成像效率。

Description

基于快速SBL算法实现分段观测ISAR高分辨成像方法及装置
技术领域
本发明属于雷达技术领域，具体涉及一种基于快速SBL算法实现分段观测ISAR高分辨成像方法及装置。
背景技术
由于逆合成孔径雷达(Inverse Synthetic Aperture Radar，ISAR)能够实现对非合作目标成像，提供目标的形状、结构、材料、运动状态和质量分布等物理信息，该独特优势使其在国防领域越来越受到世界各军事强国的重视。高分辨图像能够直观地提供更丰富的目标信息，是ISAR成像的关键技术。提高距离分辨率有两种方法，第一是通过发射具有大时宽带宽的线性调频信号，并在接收端进行匹配滤波和脉冲压缩，达到目标散射点在距离上的分辨，但这种方法成本较高。所以，一般会选择第二种方法，即通过插值或外推对稀疏子带数据进行宽带合成来提高距离分辨率，但这种方法的核心是利用阶段性缺失数据进行精确的频谱或参数估计。方位分辨率取决于相干积累角度。所谓相干积累是指脉冲之间在某一方向高度相关，使之通过一定的变换达到能量积聚，起到可分辨的效果。事实上，由于非合作目标的机动性，也很难获得较大的相干孔径数据。小视角虽简化了信号形式，提高了成像效率，但不易满足高分辨的要求。面对有限的短孔径数据，挖掘有效的先验信息和提高分辨率是ISAR成像近几年的研究热点之一。在ISAR系统中，雷达任务一般分为目标跟踪、扫描和宽带成像等模式，通过发射窄带脉冲在宽视角范围内进行波束扫描实现对目标的捕获和追踪，然后发射宽频带信号进行一定数量脉冲的相干积累以实现对目标的高分辨成像。在实际多目标成像中，单天线雷达系统需要在不同的功能之间进行不断的交替和切换，以保证雷达成像信息感知的时效性。多功能相控阵雷达的设计即为了实现雷达多项任务之间的调度。但是，对观测场景或者目标连续长时间不间断的宽带观测一般是非常困难的，目标的方位孔径观测往往出现不连续性，导致观测数据是阶段性缺失的，即稀疏孔径。在宽带组网雷达系统中，由于现有体制限制，对单一目标而言，通常每个传统器仅能获得有限时间的观测，对应目标的稀疏孔径观测。此外，在实际雷达工作中，由于外界环境干扰或自身的影响，雷达在某些特定时间内的测量观测数据是无效的。针对距离维或方位维上阶段性缺失的观测数据，或称为分段观测数据，传统的距离多普勒(RD)算法的成像结果存在很强的重影、栅瓣和副瓣问题，导致图像的分辨率下降。因此，研究分段观测ISAR高分辨成像具有重要的现实意义。
针对雷达目标回波的稀疏特性，将雷达成像模型转化为稀疏表示模型，并采用稀疏重构方法来对雷达目标参数进行优化求解，由此发展出了基于稀疏表示理论的雷达成像技术。稀疏表示理论发展至今已开发出众多的稀疏重构算法，在众多算法中，稀疏贝叶斯学习(SBL)算法具有更强的鲁棒性和更高的估计精度，故在理论和应用方面都引起了研究者的研究兴趣。SBL算法是利用信号和图像的先验信息，通过贝叶斯推理方法来求解稀疏信号参数的后验概率。在贝叶斯计算求解过程中能够自动的估计信号参数，而无需人工参数干预。而且，由于SBL算法考虑了信号的噪声统计信息，在低信噪比条件下往往可获得优于其它的稀疏重构算法的结果。
SBL算法是通过迭代实现信号重构的。但SBL算法每次迭代包含一个矩阵求逆运算，该矩阵维度与观测数据长度一样。若用传统的直接求逆方法求解，其计算复杂度与观测数据长度的立方成正比，所以，当观测数据样本数较多时，计算时间往往很长。针对这一问题，众多学者已经提出了一些快速SBL算法，但这些快速算法都采用了一些近似，会降低结果的准确性。若用于分段观测ISAR成像，结果的准确性将会更差。
发明内容
为了解决现有技术中存在的上述问题，本发明提供了一种基于快速SBL算法实现分段观测ISAR高分辨成像方法及装置。本发明要解决的技术问题通过以下技术方案实现：
本发明实施例提供了一种基于快速SBL算法实现分段观测ISAR高分辨成像方法，包括步骤：
S1、根据分段观测数据中有效采样样本和缺失采样样本对所述分段观测数据进行建模，得到分段观测数据模型；
S2、基于所述分段观测数据模型下的傅里叶字典矩阵，结合信号精度值和噪声精度值，采用快速傅里叶变换计算稀疏贝叶斯学习算法对稀疏信号进行迭代过程中的待求逆矩阵；
S3、基于所述待求逆矩阵，利用Levinson-Durbin迭代算法求解逆矩阵的G-S分解因子；
S4、结合所述逆矩阵的G-S分解因子，利用快速傅里叶算法求解所述稀疏贝叶斯学习算法对稀疏信号进行迭代过程中后验分布的均值和协方差的对角线元素；
S5、当判断所述均值不满足收敛条件时，基于所述协方差的对角线元素和所述均值计算新的信号精度值和新的噪声精度值，并返回步骤S2；当判断所述均值满足所述收敛条件时，输出所述均值。
在本发明的一个实施例中，步骤S1包括：
根据所述缺失采样样本的位置将所述分段观测数据分为若干段，则所述分段观测数据表示为：
Figure BDA0003784390080000031
所述有效数据表示为：
Figure BDA0003784390080000032
其中，
Figure BDA0003784390080000033
Ns(i)为第i段数据ys(i)的长度，Ng(i)为有效数据yg(i)的长度，Nm(i)为缺失数据ym(i)的长度，q为分段观测数据的段数。
在本发明的一个实施例中，步骤S2包括：
S21、在所述分段观测数据模型下，所述傅里叶字典矩阵中(k+1)列的傅里叶基表示为：
Figure BDA0003784390080000034
其中，ωk＝2πk/K,k＝0,...,K-1，K为超分辨倍数与观测数据的乘积；
Figure BDA0003784390080000035
为与
Figure BDA0003784390080000036
相对应的傅里叶基，
Figure BDA0003784390080000037
Figure BDA0003784390080000038
为第i段有效数据，
Figure BDA0003784390080000039
代表长度为Ng(i)的一个完整傅里叶基，
Figure BDA00037843900800000310
No(i)为基于分段观测数据模型的长度偏量；
S22、基于所述傅里叶字典矩阵，结合所述信号精度值和所述噪声精度值，采用快速傅里叶变换计算所述待求逆矩阵：
Figure BDA00037843900800000311
其中，
Figure BDA00037843900800000312
为一个维度为Ng×Ng的单位矩阵，β为噪声精度值，
Figure BDA00037843900800000313
为傅里叶字典矩阵，Λ为由1/γk按顺序构成的对角矩阵，γk为信号精度值，
Figure BDA00037843900800000314
是一个q×q的厄米特-块-托普利兹矩阵；
Figure BDA00037843900800000315
其中，Ri,j为
Figure BDA00037843900800000316
中一个维度为Ng(i)×Ng(j)的子矩阵，
Figure BDA0003784390080000041
其中，rm为Ri,j中的元素；
S23、使用两个置换矩阵和五个自定义矩阵，将所述待求逆矩阵用目标形式的待求逆矩阵表示：
Figure BDA0003784390080000042
其中，
Figure BDA0003784390080000043
和
Figure BDA0003784390080000044
为置换矩阵，
Figure BDA0003784390080000045
为自定义矩阵。
在本发明的一个实施例中，步骤S3包括：
S31、基于所述目标形式的待求逆矩阵，利用G-S分解计算所述逆矩阵：
Figure BDA0003784390080000046
其中，
Figure BDA0003784390080000047
Figure BDA0003784390080000048
Figure BDA0003784390080000049
和
Figure BDA00037843900800000410
为置换矩阵，
Figure BDA00037843900800000411
Figure BDA00037843900800000412
为自定义矩阵；
S32、计算所述逆矩阵的位移表示：
Figure BDA00037843900800000413
其中，
Figure BDA00037843900800000414
为基于
Figure BDA00037843900800000415
的自定义变量，
Figure BDA00037843900800000416
为基于
Figure BDA00037843900800000417
的自定义变量，
Figure BDA00037843900800000418
和
Figure BDA00037843900800000419
为G-S分解因子；
S33、基于所述位移表示，计算所述逆矩阵的G-S分解式：
Figure BDA00037843900800000420
其中，
Figure BDA00037843900800000421
为块托普利兹矩阵，
Figure BDA00037843900800000422
和
Figure BDA00037843900800000423
为G-S分解因子，
Figure BDA00037843900800000424
为块位移矩阵，M＝max[Ng(1),Ng(2),…,Ng(q)]，Ng(i)为有效数据yg(i)的长度；
S34、利用Levinson-Durbin迭代算法计算所述G-S分解式的所述G-S分解因子。
在本发明的一个实施例中，步骤S4包括：
S41、结合所述逆矩阵的G-S分解因子，利用快速傅里叶变换或者快速傅里叶逆变换求解所述协方差的对角线元素，所述协方差及其对角线元素构成的向量的表达式为：
Figure BDA0003784390080000051
ε＝diag(Σ)
其中，ε＝diag(Σ)代表向量ε是由矩阵Σ对角线上元素按顺序构成的，Λ为由1/γk按顺序构成的对角矩阵，γk为信号精度值，β为噪声精度值，
Figure BDA0003784390080000052
为傅里叶字典矩阵，
Figure BDA0003784390080000053
为逆矩阵；
S42、结合所述逆矩阵的G-S分解因子，利用快速傅里叶变换或者快速傅里叶逆变换求解所述均值，所述均值的表达式为：
Figure BDA0003784390080000054
其中，β为噪声精度值，
Figure BDA0003784390080000055
为傅里叶字典矩阵，
Figure BDA0003784390080000056
为逆矩阵，Λ为由1/γk按顺序构成的对角矩阵，
Figure BDA0003784390080000057
为总的有效数据。
在本发明的一个实施例中，步骤S41包括步骤：
S411、结合所述逆矩阵的G-S分解因子，利用快速傅里叶变换或者快速傅里叶逆变换求解对角线矩阵：
Figure BDA0003784390080000058
其中，
Figure BDA0003784390080000059
代表对括号中的向量做K点快速傅里叶变换，
Figure BDA00037843900800000510
和
Figure BDA00037843900800000511
为基于分段观测数据模型的自定义向量，
Figure BDA00037843900800000512
为傅里叶字典矩阵，
Figure BDA00037843900800000513
为逆矩阵；
S412、利用所述对角线矩阵、所述信号精度值和所述噪声精度值计算所述协方差：
Figure BDA00037843900800000514
其中，εk为矩阵ε的第(k+1)个值，矩阵ε是协方差矩阵对角线上的元素按顺序构成的向量，δk为对角线矩阵δ的第(k+1)个值，β为噪声精度值，γk为信号精度值。
在本发明的一个实施例中，步骤S42包括步骤：
S421、结合所述逆矩阵的G-S分解因子，利用快速傅里叶变换或者快速傅里叶逆变换求解第一向量：
Figure BDA0003784390080000061
其中，
Figure BDA0003784390080000062
为逆矩阵，
Figure BDA0003784390080000063
为总的有效数据；
S422、利用所述第一向量和所述傅里叶字典矩阵计算第二向量：
Figure BDA0003784390080000064
其中，
Figure BDA0003784390080000065
代表对括号中的向量做K点快速傅里叶逆变换，
Figure BDA0003784390080000066
为傅里叶字典矩阵，
Figure BDA0003784390080000067
为第一向量；
S423、结合所述第二向量、所述噪声精度值计算所述均值：
Figure BDA0003784390080000068
其中，β为噪声精度值，Λ为由1/γk按顺序构成的对角矩阵，
Figure BDA0003784390080000069
为第二向量。
在本发明的一个实施例中，所述收敛条件为：
Figure BDA00037843900800000610
其中，δ为设置的收敛门槛。
本发明的另一个实施例提供了一种基于快速SBL算法实现分段观测ISAR高分辨成像装置，包括：
分段观测数据模型建立模块，用于根据分段观测数据中有效采样样本和缺失采样样本对所述分段观测数据进行建模，得到分段观测数据模型；
待求逆矩阵计算模块，用于基于所述分段观测数据模型下的傅里叶字典矩阵，结合信号精度值和噪声精度值，采用快速傅里叶变换计算稀疏贝叶斯学习算法对稀疏信号进行迭代过程中的待求逆矩阵；
G-S分解因子求解模块，用于基于所述待求逆矩阵，利用Levinson-Durbin迭代算法求解逆矩阵的G-S分解因子；
协方差的对角线元素和均值求解模块，用于结合所述逆矩阵的G-S分解因子，利用快速傅里叶算法求解所述稀疏贝叶斯学习算法对稀疏信号进行迭代过程中后验分布的均值和协方差的对角线元素；
判断模块，用于当判断所述均值不满足收敛条件时，基于所述协方差的对角线元素和所述均值计算新的信号精度值和新的噪声精度值，并返回所述待求逆矩阵计算模块；当判断所述均值满足所述收敛条件时，输出所述均值。
与现有技术相比，本发明的有益效果：
本发明的成像方法采用基于傅里叶字典矩阵的稀疏贝叶斯学习算法，实现了阶段性缺失数据的高分辨率成像，成像结果准确性高；利用G-S分解表示稀疏贝叶斯学习算法迭代中的逆矩阵，避免了直接计算逆矩阵导致的高计算复杂度；涉及逆矩阵时可通过快速傅里叶变换快速求解，计算量降低了几个数量级，计算效率高；从而，该成像方法在保证成像结果准确性的基础上减小计算复杂度，大大提高了成像效率。
附图说明
图1为本发明实施例提供的一种基于快速SBL算法实现分段观测ISAR高分辨成像方法的流程示意图；
图2为本发明实施例提供的一种分段观测数据模型示意图。
图3为本发明实施例提供的一种不同算法的重构结果示意图；
图4为本发明实施例提供的一种观测数据长度不同时各种算法的性能曲线图；
图5为本发明实施例提供的一种观测数据缺失率不同时各种算法的性能曲线图；
图6为本发明实施例提供的一种观测数据所分的段数不同时各种算法的性能曲线图；
图7为本发明实施例提供的一种完整“雅克-42”数据的HRRP以及使用距离-多普勒算法和SBL算法的成像结果图；
图8为本发明实施例提供的一种分段“雅克-42”数据的HRRP以及使用距离-多普勒算法和FD-GSBL算法的成像结果图。
具体实施方式
下面结合具体实施例对本发明做进一步详细的描述，但本发明的实施方式不限于此。
实施例一
本实施例首先介绍稀疏信号重构的模型以及SBL算法。
11)稀疏信号重构是指在信号x具有一定稀疏性的前提下，根据测量数据y求解信号x的过程；其系统模型可以用一个带噪声的欠定线性系统来描述，如：
y＝Dx+e
其中，
Figure BDA0003784390080000071
是观测数据；
Figure BDA0003784390080000072
是过完备字典矩阵且K＞＞N；x指待重构的稀疏信号，即向量x中的大多数元素是零；
Figure BDA0003784390080000073
指观测噪声。
12)在众多稀疏信号重构算法中，稀疏贝叶斯学习(SBL)算法以其较高的精度引起了研究者的兴趣。这是因为SBL是一种非常重要的贝叶斯统计优化算法，它是在贝叶斯理论的基础上发展起来的，从统计的角度来实现信号重构。即在SBL框架下，待恢复信号满足一定的先验分布，然后通过贝叶斯分析得到信号的后验分布信息。SBL框架的介绍具体如下：
12a)信号的先验分布。
为了有效地提高信号的稀疏性，通常使用分层贝叶斯先验模型来描述SBL中的信号。分层模型的第一层是信号x和噪声e的建模。假设x服从零均值协方差复高斯分布Λ，e服从零均值协方差复高斯分布β-1I，则概率密度函数(PDF)为：
Figure BDA0003784390080000081
Figure BDA0003784390080000082
其中，xk代表向量x的第(k+1)个元素并且向量x中的每个元素彼此是独立的，γk代表信号精度(逆方差)，Λ是一个由1/γk按顺序构成的对角矩阵，en代表向量e的第(n+1)个元素，β代表噪声精度。
分层模型的第二层是γk和β的建模，它们都是伽马分布，其PDF分别是：
Figure BDA0003784390080000083
Figure BDA0003784390080000084
其中，Gamma(·)表示括号中的变量服从伽马分布，a和b分别是γk的形状和尺度参数，c和d分别是β的形状和尺度参数，它们被称为超参数。为获得广泛的超先验，超参数通常被设置为很小的正常数，Γ(a)表示伽马函数。
12b)信号的后验分布。
基于12a)中介绍的先验分布和测量数据y，通过使用贝叶斯公式和期望最大化(EM)算法可得到x的后验分布，其后验分布可被解析地表示为一个复高斯分布：
Figure BDA0003784390080000085
其中，该复高斯分布的协方差Σ＝(βDHD+Λ-1)-1，均值μ＝βΣDHy。
根据伍德伯里矩阵的恒等式，Σ和μ又可表示为：
Σ＝(βDHD+Λ-1)-1＝Λ-βΛDHQ-1DΛ
μ＝βΛDHQ-1y
其中，Q＝I+βDΛDH。
12c)SBL算法迭代过程。
在SBL中，信号重构是通过迭代实现的。12b)介绍的后验分布的最优均值则是重构出的信号。下面是SBL算法的迭代步骤，这里称为直接求逆的SBL(DI-SBL)：
Figure BDA0003784390080000091
Figure BDA0003784390080000092
Figure BDA0003784390080000093
Q＝I+βDΛDH
Σ(j)＝Λ-β(j)ΛDHQ-1DΛ
ε(j)＝diag(Σ(j))
μ(j)＝β(j)ΛDHQ-1y
其中，ε＝diag(Σ)代表ε是一个由矩阵Σ对角线上元素构成的向量，εk代表ε的第(k+1)个元素；
Figure BDA0003784390080000094
代表第j次迭代后得到的γk；μ和∑分别代表信号x后验概率的均值和协方差；其中μ的未知参数γk和β被称为超参数，可通过最大期望算法求解；||·||2代表l2范数。
从上述DI-SBL的迭代过程可看出SBL单次迭代过程的关键步骤是计算ε和μ，但计算过程需要求解Q-1，传统的直接求逆方法的计算复杂度与矩阵维度的立方成正比，而矩阵Q的维度与观测向量y的维度一样，若观测数据较多，计算时间往往很长，实际工程上是难以实现的。
在分段观测ISAR成像中，SBL虽可实现高分辨成像且准确性较高，但存在上述的问题，导致时间成本较高。因此，本实施例提出了一种基于傅里叶字典的快速SBL方法去实现分段观测雷达高分辨成像，该方法简称为FD-GSBL。因为在基于傅里叶字典的SBL算法中，待求逆矩阵Q是一个厄米特-块-托普利兹(Hermitian-Block-Toeplitz)矩阵，通过采用Gohberg-Semencul(G-S)分解可求解Q-1，避免了直接求逆导致的较大计算复杂度。此外，基于G-S分解因子，ε和μ可通过FFT/IFFT求解，极大地缩短了计算时间。
请参见图1，图1为本发明实施例提供的一种基于快速SBL算法实现分段观测ISAR高分辨成像方法的流程示意图。该成像方法具体包括步骤：
S1、根据分段观测数据中有效采样样本和缺失采样样本对所述分段观测数据进行建模，得到分段观测数据模型。
请参见图2，图2为本发明实施例提供的一种分段观测数据模型示意图。对于分段观测数据，其信号重构示意图如图2所示，灰框表示有效的采样样本，白框表示缺失的采样样本，可根据缺失采样样本的位置将全部分段观测数据分为q段，第i段数据ys(i)的长度为Ns(i)，有效数据yg(i)的长度为Ng(i)，缺失数据ym(i)的长度为Nm(i)。则所有分段观测数据和有效数据可分别被表示为：
Figure BDA0003784390080000101
Figure BDA0003784390080000102
其中，
Figure BDA0003784390080000103
Ns(i)为第i段数据ys(i)的长度，Ng(i)为有效数据yg(i)的长度，Nm(i)为缺失数据ym(i)的长度，q为分段观测数据的段数。
很显然，
Figure BDA0003784390080000104
和Ns(i)＝Ng(i)+Nm(i)，并且该数据的缺失率为
Figure BDA0003784390080000105
S2、基于所述分段观测数据模型下的傅里叶字典矩阵，结合信号精度值和噪声精度值，采用快速傅里叶变换计算稀疏贝叶斯学习算法对稀疏信号进行迭代过程中的待求逆矩阵。具体包括步骤：
S21、在所述分段观测数据模型下表示傅里叶字典矩阵。
首先，定义一个长度偏量No(i)：
Figure BDA0003784390080000106
因为本实施例使用的是傅里叶基构成的字典，当数据存在缺失时，该字典矩阵并不是一个完整的傅里叶字典，用
Figure BDA0003784390080000107
代表该傅里叶字典矩阵，傅里叶字典矩阵
Figure BDA0003784390080000108
中第(k+1)列的傅里叶基表示为：
Figure BDA0003784390080000109
其中，ωk＝2πk/K,k＝0,...,K-1，K为超分辨倍数与总的观测数据长度的乘积。
Figure BDA00037843900800001010
为与
Figure BDA00037843900800001011
相对应的傅里叶基。
Figure BDA0003784390080000111
表示为：
Figure BDA0003784390080000112
其中，
Figure BDA0003784390080000113
为第i段有效数据，
Figure BDA0003784390080000114
代表长度为Ng(i)的一个完整傅里叶基，No(i)为基于分段观测数据模型的长度偏量。
Figure BDA0003784390080000115
表示为：
Figure BDA0003784390080000116
S22、基于所述傅里叶字典矩阵，结合所述信号精度值和所述噪声精度值，采用快速傅里叶变换计算所述待求逆矩阵。
基于步骤S21中的傅里叶字典矩阵
Figure BDA0003784390080000117
(12b)中Q可以表示为：
Figure BDA0003784390080000118
其中，
Figure BDA0003784390080000119
为待求逆矩阵，
Figure BDA00037843900800001110
为一个维度为Ng×Ng的单位矩阵，β为噪声精度值，
Figure BDA00037843900800001111
为傅里叶字典矩阵，Λ为由1/γk按顺序构成的对角矩阵，γk为信号精度值，
Figure BDA00037843900800001112
是一个q×q的厄米特-块-托普利兹矩阵。
厄米特-块-托普利兹矩阵
Figure BDA00037843900800001113
表示为：
Figure BDA00037843900800001114
上式中，Ri,j为
Figure BDA00037843900800001115
中一个维度为Ng(i)×Ng(j)的子矩阵，其表达式为：
Figure BDA00037843900800001116
其中，rm为Ri,j中的元素，rm的计算式为：
Figure BDA00037843900800001117
从上式可看出rm可通过对1/γk做K点快速傅里叶变换(FFT)快速算得，其计算复杂度为O(Klog2K)。
基于
Figure BDA0003784390080000121
待求逆矩阵
Figure BDA0003784390080000122
可以被计算出。很明显待求逆矩阵
Figure BDA0003784390080000123
也是一个厄米特-块-托普利兹矩阵并且结构与
Figure BDA0003784390080000124
的结构一样。
S23、使用两个置换矩阵和五个自定义矩阵，将所述待求逆矩阵用目标形式的待求逆矩阵表示。
首先，将使用两个置换矩阵将待求逆矩阵
Figure BDA00037843900800001215
可被表示成两种分块矩阵的形式。
具体的，鉴于
Figure BDA0003784390080000125
中的子矩阵Qi,j是一个维度为Ng(i)×Ng(j)的托普利兹矩阵，它可以写成两种分块矩阵的形式：
Figure BDA0003784390080000126
其中，
Figure BDA0003784390080000127
Figure BDA0003784390080000128
Figure BDA0003784390080000129
代表Qi,j中一个维度为(Ng(i)-1)×(Ng(j)-1)的子矩阵，并且也是托普利兹矩阵。
基于子矩阵Qi,j的两种形式，定义下面七个矩阵，这些矩阵的元素均来自于每个Qi,j。为方便介绍，使用上标i,j代表矩阵Qi,j中的元素，例如，
Figure BDA00037843900800001210
和
Figure BDA00037843900800001211
分别代表Qi,j中的q0和q1以及
Figure BDA00037843900800001212
Figure BDA00037843900800001213
Figure BDA00037843900800001214
Figure BDA0003784390080000131
Figure BDA0003784390080000132
类似地，定义出
Figure BDA0003784390080000133
和
Figure BDA0003784390080000134
由于待求逆矩阵
Figure BDA0003784390080000135
是一个厄米特矩阵，因此可得到
Figure BDA0003784390080000136
和
Figure BDA0003784390080000137
使用上述七个自定义矩阵和两个置换矩阵，可将待求逆矩阵
Figure BDA0003784390080000138
写成下面目标形式：
Figure BDA0003784390080000139
其中，
Figure BDA00037843900800001310
和
Figure BDA00037843900800001311
为置换矩阵，
Figure BDA00037843900800001312
为自定义矩阵。
需要说明的是，上述七个自定义矩阵中，由于
Figure BDA00037843900800001313
和
Figure BDA00037843900800001314
因此，在将待求逆矩阵写目标形式时，仅使用到五个自定义矩阵。
S3、基于所述待求逆矩阵，利用Levinson-Durbin迭代算法求解逆矩阵的G-S分解因子。
S31、基于所述目标形式的待求逆矩阵，利用G-S分解计算所述逆矩阵。
具体的，对步骤S23求解出的待求逆矩阵
Figure BDA00037843900800001315
利用块矩阵求逆公式计算，得到逆矩阵：
Figure BDA00037843900800001316
Figure BDA00037843900800001317
Figure BDA00037843900800001318
Figure BDA00037843900800001319
Figure BDA00037843900800001320
其中，
Figure BDA00037843900800001321
和
Figure BDA00037843900800001322
为置换矩阵，
Figure BDA00037843900800001323
为自定义矩阵。
S32、计算所述逆矩阵的位移表示。
因为
Figure BDA0003784390080000141
是由q2个子矩阵构成的一个维度为(Ng-q)×(Ng-q)的矩阵，那
Figure BDA0003784390080000142
也可被划分成q2个子矩阵的形式，用
Figure BDA0003784390080000143
代表
Figure BDA0003784390080000144
中其中一个维度为(Ng(i)-1)×(Ng(j)-1)的子矩阵。通过使用置换矩阵
Figure BDA0003784390080000145
和
Figure BDA0003784390080000146
可构造出两个矩阵：
Figure BDA0003784390080000147
和
Figure BDA0003784390080000148
它们的其中一个子矩阵分别是：
Figure BDA0003784390080000149
和
Figure BDA00037843900800001410
接下来，定义一个维度为Ng×Ng的块位移矩阵
Figure BDA00037843900800001411
如：
Figure BDA00037843900800001412
其中，SN是一个维度为N×N的单位下三角位移矩阵，即
Figure BDA00037843900800001413
基于块位移矩阵
Figure BDA00037843900800001414
可得
Figure BDA00037843900800001415
那么，逆矩阵
Figure BDA00037843900800001416
的位移表示
Figure BDA00037843900800001417
公式表示如下：
Figure BDA00037843900800001418
令
Figure BDA0003784390080000151
Figure BDA0003784390080000152
因此，位移表示
Figure BDA0003784390080000153
为：
Figure BDA0003784390080000154
其中，
Figure BDA0003784390080000155
为基于
Figure BDA0003784390080000156
的自定义变量，
Figure BDA0003784390080000157
为基于
Figure BDA0003784390080000158
的自定义变量，
Figure BDA0003784390080000159
和
Figure BDA00037843900800001510
为G-S分解因子。从位移表示
Figure BDA00037843900800001511
的表达式可知逆矩阵
Figure BDA00037843900800001512
的位移秩为2q。
S33、基于所述位移表示，计算所述逆矩阵的G-S分解式。
具体的，使用步骤S32中的维度为Ng×Ng的块位移矩阵
Figure BDA00037843900800001513
和位移表示
Figure BDA00037843900800001514
逆矩阵
Figure BDA00037843900800001515
表示为：
Figure BDA00037843900800001516
其中，M＝max[Ng(1),Ng(2),…,Ng(q)]。
定义一个托普利兹矩阵
Figure BDA00037843900800001517
和一个块托普利兹矩阵
Figure BDA00037843900800001518
如：
Figure BDA00037843900800001519
Figure BDA00037843900800001520
其中，tN和
Figure BDA00037843900800001521
分别表示一个长度为N和Ng的向量。
那么，逆矩阵
Figure BDA00037843900800001522
的G-S分解表达式为：
Figure BDA00037843900800001523
其中，
Figure BDA00037843900800001524
为块托普利兹矩阵，
Figure BDA00037843900800001525
和
Figure BDA00037843900800001526
为G-S分解因子，
Figure BDA00037843900800001527
为块位移矩阵，M＝max[Ng(1),Ng(2),…,Ng(q)]，Ng(i)为有效数据yg(i)的长度。
上式称为逆矩阵
Figure BDA0003784390080000161
的G-S分解，
Figure BDA0003784390080000162
和
Figure BDA0003784390080000163
称为G-S分解的分解因子。
S34、利用Levinson-Durbin迭代算法计算所述G-S分解式的所述G-S分解因子。
受Levinson-Durbin(L-D)算法的影响，本实施例提出了一种迭代方法计算逆矩阵
Figure BDA0003784390080000164
的G-S分解因子。具体如下：
令
Figure BDA0003784390080000165
Figure BDA0003784390080000166
Figure BDA0003784390080000167
输入：
Figure BDA0003784390080000168
和
Figure BDA0003784390080000169
初始迭代：
Figure BDA00037843900800001610
迭代过程：
Figure BDA00037843900800001611
Figure BDA00037843900800001612
Figure BDA00037843900800001613
Figure BDA00037843900800001614
Figure BDA00037843900800001615
Figure BDA00037843900800001616
其中，
Figure BDA00037843900800001617
Figure BDA00037843900800001618
输出：
Figure BDA0003784390080000171
然后利用步骤S31中的Wq和Vq的计算式求解Wq和Vq。
S4、结合所述逆矩阵的G-S分解因子，利用快速傅里叶算法求解所述稀疏贝叶斯学习算法对稀疏信号进行迭代过程中后验分布的协方差的对角线元素和均值。
本实施例中，快速傅里叶算法包括快速傅里叶变换或者快速傅里叶逆变换。
具体的，基于傅里叶字典，协方差Σ及其对角线元素构成的向量ε和均值μ表示为：
Figure BDA0003784390080000172
ε＝diag(Σ)
Figure BDA0003784390080000173
其中，ε＝diag(Σ)代表向量ε是由矩阵Σ对角线上元素按顺序构成的，Λ为由1/γk按顺序构成的对角矩阵，γk为信号精度值，β为噪声精度值，
Figure BDA0003784390080000174
为傅里叶字典矩阵，
Figure BDA0003784390080000175
为逆矩阵；
Figure BDA0003784390080000176
为总的有效数据。
S41、结合所述逆矩阵的G-S分解因子，利用快速傅里叶变换或者快速傅里叶逆变换求解所述协方差的对角线元素。
因为在协方差公式Σ中，Λ是一个对角阵，协方差Σ的对角线元素的计算可分为两步：
Figure BDA0003784390080000177
Figure BDA0003784390080000178
其中，εk和δk分别表示对角线矩阵ε和目标矩阵δ的第(k+1)个值，
Figure BDA0003784390080000179
为傅里叶字典矩阵，
Figure BDA00037843900800001710
为逆矩阵。又因为信号精度值γk和噪声精度值β可通过(12b)计算得到，只需要快速计算目标矩阵δ。
S411、结合所述逆矩阵的G-S分解因子，利用快速傅里叶变换或者快速傅里叶逆变换求解对角线矩阵。
具体的，将步骤S21中的傅里叶字典矩阵和步骤S31中的求逆矩阵代入对角线矩阵
Figure BDA00037843900800001711
中，得：
Figure BDA00037843900800001712
Figure BDA00037843900800001713
Figure BDA0003784390080000181
其中，
Figure BDA0003784390080000182
代表
Figure BDA0003784390080000183
的第i个子向量，
Figure BDA0003784390080000184
代表
Figure BDA0003784390080000185
的第i个子向量。
Figure BDA0003784390080000186
的第(k+1)个值为：
Figure BDA0003784390080000187
其中，
Figure BDA0003784390080000188
表示矩阵
Figure BDA0003784390080000189
的第m条对角线上所有元素的和。其可通过下面的方法快速求解：
如果Ng(i)≥Ng(j)，令
Figure BDA00037843900800001810
其计算表达式为：
Figure BDA00037843900800001811
其中，若
Figure BDA00037843900800001812
则
Figure BDA00037843900800001813
注意这里的上标i指向量
Figure BDA00037843900800001814
中的元素。
很明显cij的表达式右侧是一个托普利兹矩阵与一个向量的乘积。而它们的乘积可以通过FFT/IFFT有效实现。
如果Ng(i)≥Ng(j)+2，令
Figure BDA00037843900800001815
Figure BDA00037843900800001816
的计算式为：
Figure BDA00037843900800001817
如果Ng(i)＜Ng(j)，
Figure BDA00037843900800001818
可通过相同的方法计算。
从上述
Figure BDA00037843900800001819
的计算式可看出，
Figure BDA00037843900800001820
可通过FFT快速计算：
令
Figure BDA0003784390080000191
其中，0表示使向量
Figure BDA0003784390080000192
的维度为N的一个零向量，0ij表示一个维度为
Figure BDA0003784390080000193
的零向量，0ji表示一个维度为(No(i)-No(j)-Ng(j))的零向量。然后，通过对
Figure BDA0003784390080000194
做K点的快速傅里叶变换FFT计算
Figure BDA0003784390080000195
Figure BDA0003784390080000196
可通过同样的方法求得。所以对角线矩阵δ可通过快速傅里叶变换FFT快速计算，即：
Figure BDA0003784390080000197
其中，
Figure BDA0003784390080000198
代表对括号中的向量做K点FFT(快速傅里叶变换)，
Figure BDA0003784390080000199
和
Figure BDA00037843900800001910
为基于分段观测数据模型的自定义向量；
S412、利用所述对角线矩阵、所述信号精度值和所述噪声精度值计算所述协方差的对角线元素。
具体的，在计算得到对角线矩阵δ后，通过对角线矩阵δ、噪声精度值β和1/γk的点乘计算ε，从而计算得到协方差Σ的对角线元素。即：
ε＝diag(Σ)
Figure BDA00037843900800001911
其中，εk为矩阵ε的第(k+1)个值，矩阵ε是协方差矩阵对角线上的元素按顺序构成的向量，δk为对角线矩阵δ的第(k+1)个值，β为噪声精度值，γk为信号精度值。
S42、结合所述逆矩阵的G-S分解因子，利用快速傅里叶变换或者快速傅里叶逆变换求解所述均值。具体包括步骤：
S421、结合所述G-S分解因子和所述逆矩阵，利用快速傅里叶变换或者快速傅里叶逆变换求解第一向量：
Figure BDA00037843900800001912
其中，
Figure BDA00037843900800001913
为逆矩阵，
Figure BDA00037843900800001914
为有效数据。
具体的，将步骤S31中的求逆矩阵代入
Figure BDA00037843900800001915
的表达式，可见
Figure BDA00037843900800001916
的右侧是一些托普利兹矩阵与向量的乘积，所以
Figure BDA00037843900800001917
可通过快速傅里叶变换FFT/快速傅里叶逆变换IFFT快速计算。
S422、利用所述第一向量和所述傅里叶字典矩阵计算第二向量：
Figure BDA0003784390080000201
具体的，将
Figure BDA0003784390080000202
分成q段，第i段的长度为Ng(i)。令
Figure BDA0003784390080000203
其中，
Figure BDA0003784390080000204
则
Figure BDA0003784390080000205
的计算式为：
Figure BDA0003784390080000206
其中，
Figure BDA0003784390080000207
代表对括号中的向量做K点IFFT，
Figure BDA00037843900800002013
为傅里叶字典矩阵，
Figure BDA0003784390080000208
为第一向量。
S423、结合所述第二向量、所述噪声精度值计算所述均值。
具体的，通过
Figure BDA0003784390080000209
β和1/γk的点乘计算μ：
Figure BDA00037843900800002010
其中，β为噪声精度值，Λ为由1/γk按顺序构成的对角矩阵，
Figure BDA00037843900800002011
为第二向量。
本实施例利用在基于傅里叶字典的SBL算法中，待求逆矩阵是一个Block-Toeplitz矩阵，该逆矩阵可通过G-S分解表示。同时，涉及该逆矩阵的运算可通过FFT/IFFT快速计算，在保证成像结果准确性的基础上减小计算复杂度，使得运算量按数量级程度减少，实现了较好的实时性。
S5、当判断所述均值不满足收敛条件时，基于所述协方差的对角线元素和所述均值计算新的信号精度值和新的噪声精度值，并返回步骤S2；当判断所述均值满足所述收敛条件时，输出所述均值，完成高分辨率成像。
具体的，设置收敛门槛δ，并判断每次迭代得到的μ值是否满足收敛条件：
Figure BDA00037843900800002012
若不满足收敛条件，则继续进行迭代；若满足收敛条件，则可得到最优的均值，即重构出的信号，实现高分辨成像。
本实施例通过将基于傅里叶字典的SBL算法应用到分段观测ISAR成像中，实现了阶段性缺失数据的高分辨成像，成像结果准确性高。本实例所提快速算法是用G-S分解表示SBL迭代中的逆矩阵，避免了直接计算逆矩阵导致的高计算复杂度。同时涉及该逆矩阵的式子可通过FFT快速求解，计算量降低了几个数量级，计算效率高。从而，在保证成像结果准确性的基础上减小计算复杂度，大大提高了成像效率。
实施例二
在实施例一的基础上，本实施例通过仿真实验对实施例一的基于快速SBL算法实现分段观测ISAR高分辨成像方法进行进一步说明。
SBL算法参数设置：初始值
Figure BDA0003784390080000211
β(0)＝1；超参数a＝b＝c＝d＝10-6；收敛门限δ＝10-3；频率采样因子K/Ns＝4。为了较明显地看出本发明所提出的FD-GSBL算法的成像方法的性能，在本实施例中加入现有的一些典型的稀疏信号重构方法与之进行对比，包括快速迭代自适应迭代算法(FIAA)、正交匹配追踪(OMP)、S-ESBL和DI-SBL算法。这里，S-ESBL算法是已提出的一种近似快速SBL算法，DI-SBL是指直接计算SBL算法。
模拟仿真实验：模拟观测数据来自一个有25个随机频点的模拟信号，信噪比10dB。
实测数据实验：实测观测数据来自雅克-42飞机。用于采集ISAR数据的雷达工作在c波段，频带为400mhz，脉冲重复频率为300hz。距离窗口有256个采样点，成像时间包含256个脉冲。
为了表明算法的重构性能，定义信号重构的归一化均方根误差(nRMSE)为：
Figure BDA0003784390080000212
其中，
Figure BDA0003784390080000213
代表重构出的信号值，x代表真实信号值。
1.2实验内容及结果
第一步，利用软件MATLAB R2020b对模拟观测数据进行信号重构。该模拟数据长度为512，分为4段，缺失率为50％。请参见图3，图3为本发明实施例提供的一种不同算法的重构结果示意图，其中，图3(a)是FIAA算法的重构结果图；图3(b)是OMP算法的重构结果图；图3(c)是S-ESBL算法的重构结果图；图3(d)是DI-SBL算法的重构结果图；图3(e)是FD-GSBL算法的重构结果图。
表1给出了上述各种算法信号重构的时间及归一化均方根误差。
表1各种算法信号重构的时间及归一化均方根误差对比
FIAA OMP S-ESBL DI-SBL FD-GSBL
重构时间/s 0.7910 0.0548 5.8342 28.2271 1.5619
nRMSE 0.3230 0.6772 0.4500 0.2498 0.2498
第二步，做蒙特卡洛实验，比较不同参数下算法的性能图。结果如图4、5、6所示。图4为本发明实施例提供的一种观测数据长度不同时各种算法的性能曲线图，其中图4(a)是重构的计算时间，曲线图上的时间值是取对数之后的；图4(b)图是重构的归一化均方根误差，图4(c)图是归一化均方根误差的方差。图5为本发明实施例提供的一种观测数据缺失率不同时各种算法的性能曲线图，其中图5(a)图、图5(b)图和图5(c)图分别代表计算时间、归一化均方根误差和方差。图6为本发明实施例提供的一种观测数据所分的段数不同时各种算法的性能曲线图，其中图6(a)图、图6(b)图和图6(c)图分别代表计算时间、归一化均方根误差和方差。
第三步，利用软件MATLAB R2020b对实测数据进行成像。为了与分段观测数据的成像效果做对比，本实施例提供了一种完整观测数据的成像结果图。请参见图7，图7为本发明实施例提供的一种完整“雅克-42”数据的成像结果图，其中，图7(a)是高分辨率距离像(HRRP)，图7(b)是传统距离-多普勒算法的成像结果，图7(c)是DI-SBL算法的成像结果。对于分段观测数据，假设“雅克-42”数据在方位维度上发生缺失，其MR为50％。请参见图8，图8为本发明实施例提供的一种分段“雅克-42”数据的HRRP以及使用距离-多普勒算法和FD-GSBL算法的成像结果图，其中，图8(a)为分段“雅克-42”数据的HRRP图，图8(b)为使用距离-多普勒算法的成像结果图，图8(c)为使用FD-GSBL算法的成像结果图。完整数据及缺失数据成像的距离维过采样因子为4。所有成像图的动态显示范围为40dB。
表2给出了上述分段观测实测数据成像中DI-SBL和FD-GSBL实现方式的平均运行时间。
表2分段观测实测数据成像中DI-SBL和FD-GSBL的平均运行时间对比
DI-SBL FD-GSBL
时间/s 8.9184 1.9420
1.3结果分析
从图3可以看出，当信号的两个频率值相差一个最小频率分辨率单元时，FIAA、DI-SBL和FD-GSBL算法的信号重构结果非常好，说明它们具有更高的分辨率。而S-ESBL和OMP算法的信号重构结果较差。DI-SBL和FD-GSBL算法的信号重构结果相同。此外，在表1列出的各种算法的nRMSE也证明了上述结论，即FIAA、DI-SBL和FD-GSBL算法的nRMSE相对较小，DI-SBL和FD-GSBL算法的nRMSE相同。通过比较表1中的计算时间，可以得到OMP算法的计算时间最短，FD-GSBL比DI-SBL快18倍以上。
图4、5、6展示了一些变量对算法性能的影响。从图4(a)和图4(b)可以看出，随着观测数据总长度的增加，算法的计算时间变长，nRMSE逐渐减小；对于不同的MR，数据的MR越大，有效数据越少。S-ESBL、DI-SBL和FD-GSBL算法单次迭代的时间减少，但当达到收敛阈值时，总迭代次数却增加。如图5(a)和图5(b)所示，随着MR的增加，算法的计算时间在一定范围内减小，而nRMSE增加。FD-GSBL算法的计算时间比DI-SBL算法短好几倍。此外，当MR大于40％时，S-ESBL算法的nRMSE变大，并随着MR的增加而迅速增加，这是由于S-ESBL采用了某种近似，缺失样本越多，重构效果越差；当MR大于70％时，FIAA的nRMSE变大；DI-SBL和FD-GSBL的nRMSE也有所增加，但增幅很小，即使在MR为80％时，误差值也是可以接受的。从图6(a)和图6(b)可知，q值只影响FIAA和FD-GSBL的计算复杂度，而对重构误差值没有影响。并且随着q值的增加，FIAA和FD-GSBL的计算时间变长。这是因为FIAA算法中协方差矩阵的逆矩阵以及FD-GSBL算法中所求的逆矩阵的位移秩均为2q，q的值越大，计算时间越高。此外，为比较算法的稳定性，算法的nRMSE的方差图也在图4(c)、图5(c)和图6(c)中给出。很显然，这些算法的方差值均小于0.03，说明各种算法均具有良好的稳定性。
为验证本实施例所提快速算法的有效性，用传统的距离-多普勒算法和FD-GSBL算法分别对“雅克-42”飞机的测量数据进行成像。完整测量数据的高分辨距离像及成像结果如图7所示，分段测量数据的高分辨距离像及成像结果如图8所示。从图7可以看出，距离-多普勒算法的成像结果旁瓣水平较高，DI-SBL算法的成像结果较好。从图8可以看出，与完整数据的成像结果相比，距离-多普勒算法对分段数据的成像结果具有更高的旁瓣水平，而本文提出的FD-GSBL算法的成像结果更好，说明FD-GSBL算法具有较高的成像分辨率。
表2给出了分段实测实验中DI-SBL和FD-GSBL算法的计算时间。很明显，与DI-SBL相比，FD-GSBL算法的计算时间很短。因为实验中使用的实测数据长度只有128，MR为50％，有效数据量小，所以FD-GSBL算法的加速效果不明显。总之，可看出即使在MR较大的情况下，FD-GSBL算法也能获得良好的成像结果，并且计算时间较短。
实施例三
在实施例一的基础上，本实施例还提供了一种基于快速SBL算法实现分段观测ISAR高分辨成像装置，该成像装置包括：分段观测数据模型建立模块、待求逆矩阵计算模块、G-S分解因子求解模块、协方差的对角线元素和均值求解模块、判断模块。
分段观测数据模型建立模块用于根据分段观测数据中有效采样样本和缺失采样样本对所述分段观测数据进行建模，得到分段观测数据模型。待求逆矩阵计算模块用于基于所述分段观测数据模型下的傅里叶字典矩阵，结合信号精度值和噪声精度值，采用快速傅里叶变换计算稀疏贝叶斯学习算法对稀疏信号进行迭代过程中的待求逆矩阵。G-S分解因子求解模块用于基于所述待求逆矩阵，利用Levinson-Durbin迭代算法求解逆矩阵的G-S分解因子。协方差的对角线元素和均值求解模块用于结合所述逆矩阵的G-S分解因子，利用快速傅里叶算法求解所述稀疏贝叶斯学习算法对稀疏信号进行迭代过程中后验分布的协方差的对角线元素和均值。判断模块用于当判断所述均值不满足收敛条件时，基于所述协方差的对角线元素和所述均值计算新的信号精度值和新的噪声精度值，并返回所述待求逆矩阵计算模块；当判断所述均值满足所述收敛条件时，输出所述均值。
本发明实施例提供的基于快速SBL算法实现分段观测ISAR高分辨成像装置，可以执行上述方法实施例，其实现原理和技术效果类似，在此不再赘述。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (9)

1.一种基于快速SBL算法实现分段观测ISAR高分辨成像方法，其特征在于，包括步骤：
S1、根据分段观测数据中有效采样样本和缺失采样样本对所述分段观测数据进行建模，得到分段观测数据模型；
S2、基于所述分段观测数据模型下的傅里叶字典矩阵，结合信号精度值和噪声精度值，采用快速傅里叶变换计算稀疏贝叶斯学习算法对稀疏信号进行迭代过程中的待求逆矩阵；
S3、基于所述待求逆矩阵，利用Levinson-Durbin迭代算法求解逆矩阵的G-S分解因子；
S4、结合所述逆矩阵的G-S分解因子，利用快速傅里叶算法求解所述稀疏贝叶斯学习算法对稀疏信号进行迭代过程中后验分布的均值和协方差的对角线元素；
S5、当判断所述均值不满足收敛条件时，基于所述协方差的对角线元素和所述均值计算新的信号精度值和新的噪声精度值，并返回步骤S2；当判断所述均值满足所述收敛条件时，输出所述均值。
2.根据权利要求1所述的基于快速SBL算法实现分段观测ISAR高分辨成像方法，其特征在于，步骤S1包括：
根据所述缺失采样样本的位置将所述分段观测数据分为若干段，则所述分段观测数据表示为：
Figure FDA0003784390070000011
有效数据表示为：
Figure FDA0003784390070000012
其中，
Figure FDA0003784390070000013
Ns(i)为第i段数据ys(i)的长度，Ng(i)为有效数据yg(i)的长度，Nm(i)为缺失数据ym(i)的长度，q为分段观测数据的段数。
3.根据权利要求1所述的基于快速SBL算法实现分段观测ISAR高分辨成像方法，其特征在于，步骤S2包括：
S21、在所述分段观测数据模型下，所述傅里叶字典矩阵中(k+1)列的傅里叶基表示为：
Figure FDA0003784390070000021
其中，ωk＝2πk/K,k＝0,...,K-1，K为超分辨倍数与分段观测数据总长度的乘积；
Figure FDA0003784390070000022
为与
Figure FDA0003784390070000023
相对应的傅里叶基，
Figure FDA0003784390070000024
Figure FDA0003784390070000025
为第i段有效数据，
Figure FDA0003784390070000026
代表长度为Ng(i)的一个完整傅里叶基，
Figure FDA0003784390070000027
No(i)为基于分段观测数据模型的长度偏量；
S22、基于所述傅里叶字典矩阵，结合所述信号精度值和所述噪声精度值，采用快速傅里叶变换计算所述待求逆矩阵：
Figure FDA0003784390070000028
其中，
Figure FDA0003784390070000029
为维度为Ng×Ng的单位矩阵，β为噪声精度值，
Figure FDA00037843900700000210
为傅里叶字典矩阵，Λ为由1/γk按顺序构成的对角矩阵，γk为信号精度值，
Figure FDA00037843900700000211
是一个q×q的厄米特-块-托普利兹矩阵；
Figure FDA00037843900700000212
其中，Ri,j为
Figure FDA00037843900700000213
中一个维度为Ng(i)×Ng(j)的子矩阵，
Figure FDA00037843900700000214
其中，rm为Ri,j中的元素；
S23、使用两个置换矩阵和五个自定义矩阵，将所述待求逆矩阵用目标形式的待求逆矩阵表示：
Figure FDA0003784390070000031
其中，
Figure FDA0003784390070000032
和
Figure FDA0003784390070000033
为置换矩阵，
Figure FDA0003784390070000034
为自定义矩阵。
4.根据权利要求3所述的基于快速SBL算法实现分段观测ISAR高分辨成像方法，其特征在于，步骤S3包括：
S31、基于所述目标形式的待求逆矩阵，利用G-S分解计算所述逆矩阵：
Figure FDA0003784390070000035
其中
Figure FDA0003784390070000036
Figure FDA0003784390070000037
Figure FDA0003784390070000038
和
Figure FDA0003784390070000039
为置换矩阵，
Figure FDA00037843900700000310
Figure FDA00037843900700000311
为自定义矩阵；
S32、计算所述逆矩阵的位移表示：
Figure FDA00037843900700000312
其中，
Figure FDA00037843900700000313
为基于
Figure FDA00037843900700000314
的自定义变量，
Figure FDA00037843900700000315
为基于
Figure FDA00037843900700000316
的自定义变量，
Figure FDA00037843900700000317
和
Figure FDA00037843900700000318
为G-S分解因子；
S33、基于所述位移表示，计算所述逆矩阵的G-S分解式：
Figure FDA00037843900700000319
其中，
Figure FDA00037843900700000320
为块托普利兹矩阵，
Figure FDA00037843900700000321
和
Figure FDA00037843900700000322
为G-S分解因子，
Figure FDA00037843900700000323
为块位移矩阵，M＝max[Ng(1),Ng(2),…,Ng(q)]，Ng(i)为有效数据yg(i)的长度；
S34、利用Levinson-Durbin迭代算法计算所述G-S分解式的所述G-S分解因子。
5.根据权利要求1所述的基于快速SBL算法实现分段观测ISAR高分辨成像方法，其特征在于，步骤S4包括：
S41、结合所述逆矩阵的G-S分解因子，利用快速傅里叶变换或者快速傅里叶逆变换求解所述协方差的对角线元素，所述协方差及其对角线元素构成的向量表达式为：
Figure FDA0003784390070000041
ε＝diag(Σ)
其中，ε＝diag(Σ)代表向量ε是由矩阵Σ对角线上元素按顺序构成的，Λ为由1/γk按顺序构成的对角矩阵，γk为信号精度值，β为噪声精度值，
Figure FDA0003784390070000042
为傅里叶字典矩阵，
Figure FDA0003784390070000043
为逆矩阵；
S42、结合所述逆矩阵的G-S分解因子，利用快速傅里叶变换或者快速傅里叶逆变换求解所述均值，所述均值的表达式为：
Figure FDA0003784390070000044
其中，β为噪声精度值，
Figure FDA0003784390070000045
为傅里叶字典矩阵，
Figure FDA0003784390070000046
为逆矩阵，Λ为由1/γk按顺序构成的对角矩阵，
Figure FDA0003784390070000047
为总的有效数据。
6.根据权利要求5所述的基于快速SBL算法实现分段观测ISAR高分辨成像方法，其特征在于，步骤S41包括步骤：
S411、结合所述逆矩阵的G-S分解因子，利用快速傅里叶变换或者快速傅里叶逆变换求解对角线矩阵：
Figure FDA0003784390070000048
其中，
Figure FDA0003784390070000049
代表对括号中的向量做K点快速傅里叶变换，
Figure FDA00037843900700000410
和
Figure FDA00037843900700000411
为基于分段观测数据模型的自定义向量，
Figure FDA00037843900700000412
为傅里叶字典矩阵，
Figure FDA00037843900700000413
为逆矩阵；
S412、利用所述对角线矩阵、所述信号精度值和所述噪声精度值计算所述协方差的对角线元素：
Figure FDA0003784390070000051
其中，εk为矩阵ε的第(k+1)个值，矩阵ε是协方差矩阵对角线上的元素按顺序构成的向量，δk为对角线矩阵δ的第(k+1)个值，β为噪声精度值，γk为信号精度值。
7.根据权利要求5所述的基于快速SBL算法实现分段观测ISAR高分辨成像方法，其特征在于，步骤S42包括步骤：
S421、结合所述逆矩阵的G-S分解因子，利用快速傅里叶变换或者快速傅里叶逆变换求解第一向量：
Figure FDA0003784390070000052
其中，
Figure FDA0003784390070000053
为逆矩阵，
Figure FDA0003784390070000054
为有效数据；
S422、利用所述第一向量和所述傅里叶字典矩阵计算第二向量：
Figure FDA0003784390070000055
其中，
Figure FDA00037843900700000512
表示对括号中的向量做K点快速傅里叶逆变换，
Figure FDA0003784390070000057
为傅里叶字典矩阵，
Figure FDA0003784390070000058
为第一向量；
S423、结合所述第二向量、所述噪声精度值计算所述均值：
Figure FDA0003784390070000059
其中，β为噪声精度值，Λ为由1/γk按顺序构成的对角矩阵，
Figure FDA00037843900700000510
为第二向量。
8.根据权利要求1所述的基于快速SBL算法实现分段观测ISAR高分辨成像方法，其特征在于，所述收敛条件为：
Figure FDA00037843900700000511
其中，δ为设置的收敛门槛。
9.一种基于快速SBL算法实现分段观测ISAR高分辨成像装置，其特征在于，包括：
分段观测数据模型建立模块，用于根据分段观测数据中有效采样样本和缺失采样样本对所述分段观测数据进行建模，得到分段观测数据模型；
待求逆矩阵计算模块，用于基于所述分段观测数据模型下的傅里叶字典矩阵，结合信号精度值和噪声精度值，采用快速傅里叶变换计算稀疏贝叶斯学习算法对稀疏信号进行迭代过程中的待求逆矩阵；
G-S分解因子求解模块，用于基于所述待求逆矩阵，利用Levinson-Durbin迭代算法求解逆矩阵的G-S分解因子；
协方差的对角线元素和均值求解模块，用于结合所述逆矩阵的G-S分解因子，利用快速傅里叶算法求解所述稀疏贝叶斯学习算法对稀疏信号进行迭代过程中后验分布的均值和协方差的对角线元素；
判断模块，用于当判断所述均值不满足收敛条件时，基于所述协方差的对角线元素和所述均值计算新的信号精度值和新的噪声精度值，并返回所述待求逆矩阵计算模块；当判断所述均值满足所述收敛条件时，输出所述均值。
CN202210944216.8A 2022-08-05 2022-08-05 基于快速sbl算法实现分段观测isar高分辨成像方法及装置 Pending CN115453528A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210944216.8A CN115453528A (zh) 2022-08-05 2022-08-05 基于快速sbl算法实现分段观测isar高分辨成像方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210944216.8A CN115453528A (zh) 2022-08-05 2022-08-05 基于快速sbl算法实现分段观测isar高分辨成像方法及装置

Publications (1)

Publication Number Publication Date
CN115453528A true CN115453528A (zh) 2022-12-09

Family

ID=84296098

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210944216.8A Pending CN115453528A (zh) 2022-08-05 2022-08-05 基于快速sbl算法实现分段观测isar高分辨成像方法及装置

Country Status (1)

Country Link
CN (1) CN115453528A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116540203A (zh) * 2023-07-04 2023-08-04 西安电子科技大学 基于快速稀疏贝叶斯的宽带雷达高速目标的相参积累方法
CN116626646A (zh) * 2023-07-21 2023-08-22 西安电子科技大学 基于时频非均采样的雷达目标无网格化损失相参积累方法
CN118191756A (zh) * 2024-05-20 2024-06-14 中国空气动力研究与发展中心计算空气动力研究所 一种雷达信号检测方法及系统

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108008385A (zh) * 2017-11-20 2018-05-08 西安电子科技大学 基于稀疏贝叶斯学习的干扰环境isar高分辨成像方法
CN113030972A (zh) * 2021-04-26 2021-06-25 西安电子科技大学 基于快速稀疏贝叶斯学习的机动目标isar成像方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108008385A (zh) * 2017-11-20 2018-05-08 西安电子科技大学 基于稀疏贝叶斯学习的干扰环境isar高分辨成像方法
CN113030972A (zh) * 2021-04-26 2021-06-25 西安电子科技大学 基于快速稀疏贝叶斯学习的机动目标isar成像方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
HONGCHAO LIU 等: "Superresolution ISAR Imaging Based on Sparse Bayesian Learning", IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, 6 December 2013 (2013-12-06) *
S.M. KAY 等: "Spectrum analysis—A modern perspective", PROCEEDINGS OF THE IEEE, 30 November 1981 (1981-11-30) *
SHUANGHUI ZHANG 等: "Bayesian Bistatic ISAR Imaging for Targets With Complex Motion Under Low SNR Condition", IEEE TRANSACTIONS ON IMAGE PROCESSING, 7 February 2018 (2018-02-07) *
XIAOXIU ZHU 等: "Scence Segmentation of Multi-Band ISAR Fusion Imaging Based on MB-PCSBL", IEEE SENSORS JOURNAL, 1 February 2021 (2021-02-01) *
刘红超: "高分辨ISAR成像新方法研究", 中国优秀博士论文全文数据库 信息科技辑, 15 December 2014 (2014-12-15) *
汪玲;朱兆达;朱岱寅;: "机载ISAR舰船侧视和俯视成像时间段选择", 电子与信息学报, no. 12, 15 December 2008 (2008-12-15) *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116540203A (zh) * 2023-07-04 2023-08-04 西安电子科技大学 基于快速稀疏贝叶斯的宽带雷达高速目标的相参积累方法
CN116540203B (zh) * 2023-07-04 2023-09-22 西安电子科技大学 基于快速稀疏贝叶斯的宽带雷达高速目标的相参积累方法
CN116626646A (zh) * 2023-07-21 2023-08-22 西安电子科技大学 基于时频非均采样的雷达目标无网格化损失相参积累方法
CN116626646B (zh) * 2023-07-21 2023-09-22 西安电子科技大学 基于时频非均采样的雷达目标无网格化损失相参积累方法
CN118191756A (zh) * 2024-05-20 2024-06-14 中国空气动力研究与发展中心计算空气动力研究所 一种雷达信号检测方法及系统

Similar Documents

Publication Publication Date Title
CN109212526B (zh) 用于高频地波雷达的分布式阵列目标角度测量方法
CN110208735B (zh) 一种基于稀疏贝叶斯学习的相干信号doa估计方法
CN104749553B (zh) 基于快速稀疏贝叶斯学习的波达方向角估计方法
CN115453528A (zh) 基于快速sbl算法实现分段观测isar高分辨成像方法及装置
CN110244303B (zh) 基于sbl-admm的稀疏孔径isar成像方法
CN109298383B (zh) 一种基于变分贝叶斯推断的互质阵波达方向角估计方法
CN107271993B (zh) 一种基于最大后验的扫描雷达角超分辨成像方法
CN109116293B (zh) 一种基于离格稀疏贝叶斯的波达方向估计方法
CN108896990B (zh) 一种利用耦合模式字典学习的建筑物墙体成像方法和装置
CN109507666B (zh) 基于离网变分贝叶斯算法的isar稀疏频带成像方法
CN107621635B (zh) 一种前视海面目标角超分辨方法
CN112612006B (zh) 基于深度学习的机载雷达非均匀杂波抑制方法
CN111046591B (zh) 传感器幅相误差与目标到达角度的联合估计方法
CN111856465B (zh) 一种基于稀疏约束的前视海面目标角超分辨方法
CN110726992A (zh) 基于结构稀疏和熵联合约束的sa-isar自聚焦法
CN109298382A (zh) 一种基于期望极大算法的非均匀直线阵波达方向角估计方法
CN105699950A (zh) 基于自适应迭代前后向平滑共轭梯度的雷达杂波抑制方法
CN108646247A (zh) 基于伽马过程线性回归的逆合成孔径雷达成像方法
CN112766304A (zh) 一种基于稀疏贝叶斯学习的机动阵列方位估计方法
CN115963494A (zh) 基于快速sbl算法的周期性分段观测isar高分辨成像方法
CN112147608A (zh) 一种快速高斯网格化非均匀fft穿墙成像雷达bp方法
CN105929397A (zh) 基于lq 正则化的偏置相位中心天线成像方法
CN108919263B (zh) 基于最大互信息准则的isar高分辨成像方法
CN114609594A (zh) 非均匀性杂波中知识辅助宽带雷达目标检测器及设计方法
CN116577749A (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