CN109932717B - 基于环境统计建模的isar高分辨成像方法 - Google Patents

Abstract

本发明提出了一种基于环境统计建模的ISAR高分辨成像方法，用于解决现有技术中当环境存在非高斯强干扰及噪声时无法实现ISAR聚焦成像的技术问题。实现步骤为：获取运动目标的缺损回波矩阵S_r；获取有效回波矩阵S_e；获取实转置回波矩阵S；构造实傅里叶字典Φ；构建缺损回波矩阵S_r每个距离单元的稀疏信号表征模型；构建干扰向量ε的高斯混合模型；构建权向量ω_q的伽马过程‑高斯层级先验；利用最大后验‑期望最大MAP‑EM算法计算S_r每个距离单元散射点分布的权向量；构造复权值矩阵W_c；获取运动目标的二维逆合成孔径雷达ISAR成像结果。本发明能够在非高斯强干扰及噪声环境下获得聚焦良好的ISAR像，可用于空间与空中目标监视。

Description

基于环境统计建模的ISAR高分辨成像方法

技术领域

本发明属于雷达信号处理技术领域，涉及一种逆合成孔径雷达ISAR高分辨成像方法，具体涉及一种基于环境统计建模的ISAR高分辨成像方法，可用于非高斯干扰及噪声条件下的空间与空中目标监视。

背景技术

由于具有全天时、全天候和高分辨率等特点，逆合成孔径雷达ISAR目前已在军事和民用领域获得了广泛应用。逆合成孔径雷达ISAR通过发射宽带信号获得高距离分辨率，利用雷达与目标间的相对运动获得高方位分辨率，进而获得目标的二维高分辨图像。传统逆合成孔径雷达ISAR成像算法主要基于傅里叶变换，对应分辨率较低，旁瓣较高，并且不适用于缺损回波成像。随后，基于现代谱估计的方法使成像分辨率有了显著提高，但仅适用于高信噪比条件。

为了在低信噪比、回波缺损条件下实现逆合成孔径雷达ISAR高分辨成像，近年来又提出基于稀疏信号重构理论的ISAR高分辨成像方法。例如，李文静，陈红卫在其发表的论文“一种基于压缩感知的ISAR成像方法”(《计算机仿真》，2015，32(8)：10-13，62)中，公开了一种基于压缩感知的ISAR成像方法，该方法首先对回波信号进行去斜处理，通过二维稀疏采样得到观测信号，进而采用正交匹配追踪OMP算法对观测信号重构，最终实现逆合成孔径雷达ISAR成像。虽然该方法能够在回波缺损的情况下实现二维高分辨成像，但在较低信噪比条件下，成像质量较差。

又如，王天云，陆新飞，孙麟，陈畅，陈卫东在其发表的论文“基于贝叶斯压缩感知的ISAR自聚焦成像”(《电子与信息学报》，2015，(11)：2719-2726)中，公开了一种基于稀疏贝叶斯学习的ISAR自聚焦成像方法，该方法首先对散射点分布建立伽马-高斯层级先验，进而依据最大后验MAP准则求解模型参数，最终实现逆合成孔径雷达ISAR成像。虽然该方法能够在较低信噪比条件下实现缺损回波的二维成像，但由于假设回波中仅存在高斯噪声，因此当观测环境中存在非高斯的强干扰时，很难获得聚焦良好的逆合成孔径雷达ISAR像。

发明内容

本发明的目的在于克服上述现有技术中存在的缺陷，提出一种基于环境统计建模的ISAR高分辨成像方法，用于解决当观测环境中存在非高斯强干扰时，现有技术无法获得聚焦良好的逆合成孔径雷达ISAR像的问题。

本发明的技术思路是：根据回波信号模型，首先将逆合成孔径雷达ISAR成像问题转化为稀疏信号表征问题，接着对非高斯的干扰向量建立高斯混合模型，并对每个距离单元散射点分布的权向量建立伽马过程-高斯层级先验，进而采用最大后验-期望最大MAP-EM算法求解每个距离单元散射点分布的权向量，最终实现强干扰环境下的高分辨二维成像，具体实现步骤为：

(1)获取运动目标的缺损回波矩阵S_r：

逆合成孔径雷达ISAR向运动目标发射线性调频信号，并接收干扰环境下运动目标维数为N_r×N_a的缺损回波矩阵S_r，其中，N_r表示缺损回波矩阵S_r的距离单元数，N_a表示缺损回波矩阵S_r的方位单元数；

(2)获取有效回波矩阵S_e：

沿距离向对缺损回波矩阵S_r进行脉冲压缩，得到距离脉冲压缩后的缺损回波矩阵S_pc，并剔除S_pc中缺损的列向量，得到维数为N_r×N的有效回波矩阵S_e，其中，N表示有效方位单元数；

(3)获取实转置回波矩阵S：

对有效回波矩阵S_e进行转置，得到复转置回波矩阵S_c，并通过S_c构造维数为2N×N_r的实转置回波矩阵S：

其中，Re(·)表示取实部操作，Im(·)表示取虚部操作；

(4)构造实傅里叶字典Φ：

(4a)以

为元素，构造维数为N_a×N_a的复傅里叶字典Φ_c，并剔除Φ_c中与缺损回波矩阵S_r的缺损列序号相对应的行，得到维数为N×N_a的有效傅里叶字典Φ_e，其中，e^(·)表示以自然常数为底的指数操作，j表示虚数单位，π表示圆周率，m表示复傅里叶字典Φ_c行的序号，n表示复傅里叶字典Φ_c列的序号，行序号m与列序号n的取值范围均为[-N_a/2，N_a/2-1]；

(4b)通过Φ_e构造实傅里叶字典Φ：

(5)构建缺损回波矩阵S_r每个距离单元的稀疏信号表征模型：

通过实傅里叶字典Φ和实转置回波矩阵S，构建缺损回波矩阵S_r中第q个距离单元的稀疏信号表征模型：

S_f＝Φω_q+ε

其中，S_f表示实转置回波矩阵S的第f个列向量，ω_q表示缺损回波矩阵S_r第q个距离单元散射点分布的权向量，ε表示干扰向量，f＝q＝1，…，N_r；

(6)构建干扰向量ε的高斯混合模型：

构建干扰向量ε的高斯混合模型，该模型包含干扰向量ε的第i个元素ε_i的概率分布p(ε_i|π，μ，α)、精度向量α的第j个元素α_j的先验分布p(α_j)和混合系数π的先验分布p(π)，其中：

干扰向量ε的第i个元素ε_i的概率分布p(ε_i|π，μ，α)的表达式为：

其中，J表示干扰向量ε的聚类数，π_j表示混合系数π的第j个元素，μ_j表示均值向量μ的第j个元素，α_j表示精度向量α的第j个元素，N(·)表示高斯分布的概率密度，∑表示求和操作，i＝1，…，2N，j＝1，…，J；

精度向量α的第j个元素α_j的先验分布p(α_j)的表达式为：

p(α_j)＝Gam(α_j|a，b)

其中，a和b表示精度向量α的第j个元素α_j的超参数，Gam(·)表示伽马分布的概率密度；

混合系数π的先验分布p(π)的表达式为：

p(π)＝Dir(π|γ₁，...，γ_J)

其中，γ₁，...，γ_J表示混合系数π的超参数，Dir(·)表示狄利克雷分布的概率密度；

(7)构建权向量ω_q的伽马过程-高斯层级先验：

构建权向量ω_q的伽马过程-高斯层级先验，该先验包含权向量ω_q的先验分布p(ω_q)和协方差矩阵∑第d个对角线元素η_d的先验分布p(η_d)，其中：

权向量ω_q的先验分布p(ω_q)的表达式为：

p(ω_q)＝N(ω_q|0，∑)

其中，协方差矩阵∑为对角阵；

协方差矩阵∑第d个对角线元素η_d的先验分布p(η_d)的表达式为：

p(η_d)＝Gam(η_d|P，1)

其中，P表示协方差矩阵∑第d个对角线元素η_d的超参数，d＝1，…，2N_a；

(8)利用最大后验-期望最大MAP-EM算法计算S_r每个距离单元散射点分布的权向量：

(8a)令S_r初始距离单元序号为q＝1；

(8b)令初始迭代次数为t＝1，并初始化参数：

在S_r第q个距离单元中，令干扰向量ε的聚类数J的初始值为10，令均值向量μ包含元素的个数为10，且各元素的初始值为随机值，令精度向量α包含元素的个数为10，且各元素的初始值均为0.01，令协方差矩阵∑的初始值为对角线元素均为100的维数为2N_a×2N_a的对角阵，令混合系数π的超参数γ₁，...，γ_J均为1，令精度向量α每个元素的超参数a和b均为10^-4，令协方差矩阵∑每个对角线元素的超参数P均为1/2N_a，令类别矩阵L的初始值的维数为2N×10，且类别矩阵L每行中仅有一个元素为1，其余元素均为0，每行非零元素的列序号在整数[1，10]中随机产生；

(8c)计算干扰向量ε第j个聚类的有效数目

混合系数π^(t)、S_r第q个距离单元散射点分布的权向量

均值向量μ^(t)、精度向量α^(t)、类别矩阵L^(t)和协方差期望

的值；

(8d)更新类别矩阵L^(t)、干扰向量ε的聚类数J^(t)、混合系数π^(t)、均值向量μ^(t)和精度向量α^(t)；

(8e)判断当前迭代次数t＞1是否成立，若是，执行步骤(8f)；否则，令t＝t+1，并执行步骤(8c)；

(8f)判断

是否成立，若是，则S_r第q个距离单元散射点分布的权向量

并执行步骤(8g)；否则，令t＝t+1，并执行步骤(8c)，其中，||·||₂表示2-范数操作；

(8g)判断S_r当前距离单元序号q＝N_r是否成立，若是，得到S_r每个距离单元散射点分布的权向量

否则，令q＝q+1，并执行步骤(8b)；

(9)构造复权值矩阵W_c：

将S_r每个距离单元散射点分布的权向量

按行拼接，得到维数为2N_a×N_r的实权值矩阵W，并通过W构造复权值矩阵W_c：

W_c＝W₁+jW₂

其中，W₁表示实权值矩阵W第1，…，N_a行元素构成的矩阵，W₂表示实权值矩阵W第N_a+1，…，2N_a行元素构成的矩阵；

(10)获取运动目标的二维逆合成孔径雷达ISAR成像结果：

对复权值矩阵W_c进行转置，得到干扰环境下运动目标的二维高分辨逆合成孔径雷达ISAR成像结果。

本发明与现有技术相比，具有以下优点：

(1)本发明对干扰向量建立高斯混合模型，能够对任意形式的统计分布进行逼近，克服了现有技术对干扰向量所建模型不够准确，成像效果差的问题；

(2)本发明对每个距离单元散射点分布的权向量建立伽马过程-高斯层级先验，克服了现有技术对权向量所建模型稀疏表征能力不足、对数据表示不够灵活的问题，在干扰形式未知，低信干(噪)比环境下能获得聚焦良好的ISAR像。

附图说明

图1为本发明的实现流程图；

图2为本发明的仿真图，其中图2(a)为距离脉冲压缩后的缺损回波的仿真图，图2(b)为利用现有技术中的正交匹配追踪OMP方法对缺损回波重构后的成像结果仿真图，图2(c)为利用现有技术中的相关向量机RVM对缺损回波重构后的成像结果仿真图，图2(d)为利用本发明对缺损回波重构后的成像结果仿真图。

具体实施方式

下面结合附图和具体实施例，对本发明作进一步详细描述。

参照图1，一种基于环境统计建模的ISAR高分辨成像方法，包括如下步骤：

步骤1)获取运动目标的缺损回波矩阵S_r：

逆合成孔径雷达ISAR向运动目标发射线性调频信号，并接收干扰环境下运动目标维数为N_r×N_a的缺损回波矩阵S_r，其中，N_r表示缺损回波矩阵S_r的距离单元数，N_a表示缺损回波矩阵S_r的方位单元数，在本实施例中，N_r＝256，N_a＝512；

步骤2)获取有效回波矩阵S_e：

沿距离向对缺损回波矩阵S_r进行脉冲压缩，得到距离脉冲压缩后的缺损回波矩阵S_pc，并剔除S_pc中缺损的列向量，得到维数为N_r×N的有效回波矩阵S_e，其中，N表示有效方位单元数，在本实施例中，N＝256；

其中，距离脉冲压缩后的缺损回波矩阵S_pc，实现步骤为：

(2a)将逆合成孔径雷达ISAR全场景中心的距离作为参考距离，并选取载频和调频率与逆合成孔径雷达ISAR发射信号相同，且距离为参考距离的线性调频信号作为参考信号S_ref；

(2b)将参考信号S_ref取共轭后与接收的缺损回波矩阵S_r相乘，得到解线频调后的信号S_de，并对S_de沿快时间维作傅里叶变换，得到距离脉冲压缩后的缺损回波矩阵S_pc；

步骤3)获取实转置回波矩阵S：

对有效回波矩阵S_e进行转置，得到复转置回波矩阵S_c，并通过S_c构造维数为2N×N_r的实转置回波矩阵S：

其中，Re(·)表示取实部操作，Im(·)表示取虚部操作；

步骤4)构造实傅里叶字典Φ：

(4a)以

为元素，构造维数为N_a×N_a的复傅里叶字典Φ_c，并剔除Φ_c中与缺损回波矩阵S_r的缺损列序号相对应的行，得到维数为N×N_a的有效傅里叶字典Φ_e，其中，e^(·)表示以自然常数为底的指数操作，j表示虚数单位，π表示圆周率，m表示复傅里叶字典Φ_c行的序号，n表示复傅里叶字典Φ_c列的序号，行序号m与列序号n的取值范围均为[-N_a/2，N_a/2-1]；

(4b)通过Φ_e构造实傅里叶字典Φ：

步骤5)构建缺损回波矩阵S_r每个距离单元的稀疏信号表征模型：

通过实傅里叶字典Φ和实转置回波矩阵S，构建缺损回波矩阵S_r中第q个距离单元的稀疏信号表征模型：

S_f＝Φω_q+ε

其中，S_f表示实转置回波矩阵S的第f个列向量，ω_q表示缺损回波矩阵S_r第q个距离单元散射点分布的权向量，ε表示干扰向量，f＝q＝1，…，N_r；

步骤6)构建干扰向量ε的高斯混合模型：

由于干扰向量ε的统计特性未知，而高斯混合模型能够灵活地近似各种概率分布，因此构建干扰向量ε的高斯混合模型，该模型包含干扰向量ε的第i个元素ε_i的概率分布p(ε_i|π，μ，α)、精度向量α的第j个元素α_j的先验分布p(α_j)和混合系数π的先验分布p(π)，其中：

干扰向量ε的第i个元素ε_i的概率分布p(ε_i|π，μ，α)的表达式为：

其中，J表示干扰向量ε的聚类数，π_j表示混合系数π的第j个元素，μ_j表示均值向量μ的第j个元素，α_j表示精度向量α的第j个元素，N(·)表示高斯分布的概率密度，∑表示求和操作，i＝1，…，2N，j＝1，…，J；

精度向量α的第j个元素α_j的先验分布p(α_j)的表达式为：

p(α_j)＝Gam(α_j|a，b)

其中，a和b表示精度向量α的第j个元素α_j的超参数，Gam(·)表示伽马分布的概率密度；

混合系数π的先验分布p(π)的表达式为：

p(π)＝Dir(π|γ₁，...，γ_J)

其中，γ₁，...，γ_J表示混合系数π的超参数，Dir(·)表示狄利克雷分布的概率密度；

步骤7)构建权向量ω_q的伽马过程-高斯层级先验：

为了使权向量ω_q更加稀疏、模型更加灵活，构建权向量ω_q的伽马过程-高斯层级先验，该先验包含权向量ω_q的先验分布p(ω_q)和协方差矩阵∑第d个对角线元素η_d的先验分布p(η_d)，其中：

权向量ω_q的先验分布p(ω_q)的表达式为：

p(ω_q)＝N(ω_q|0，∑)

其中，协方差矩阵∑为对角阵；

协方差矩阵∑第d个对角线元素η_d的先验分布p(η_d)的表达式为：

p(η_d)＝Gam(η_d|P，1)

其中，P表示协方差矩阵∑第d个对角线元素η_d的超参数，d＝1，…，2N_a；

步骤8)利用最大后验-期望最大MAP-EM算法计算S_r每个距离单元散射点分布的权向量：

(8a)令S_r初始距离单元序号为q＝1；

(8b)令初始迭代次数为t＝1，并初始化参数：

在S_r第q个距离单元中，令干扰向量ε的聚类数J的初始值为10，令均值向量μ包含元素的个数为10，且各元素的初始值为随机值，令精度向量α包含元素的个数为10，且各元素的初始值均为0.01，令协方差矩阵∑的初始值为对角线元素均为100的维数为2N_a×2N_a的对角阵，令混合系数π的超参数γ₁，...，γ_J均为1，令精度向量α每个元素的超参数a和b均为10^-4，令协方差矩阵∑每个对角线元素的超参数P均为1/2N_a，令类别矩阵L的初始值的维数为2N×10，且类别矩阵L每行中仅有一个元素为1，其余元素均为0，每行非零元素的列序号在整数[1，10]中随机产生；

(8c)计算干扰向量ε第j个聚类的有效数目

混合系数π^(t)、S_r第q个距离单元散射点分布的权向量

均值向量μ^(t)、精度向量α^(t)、类别矩阵L^(t)和协方差期望

的值，计算公式分别为：

干扰向量ε第j个聚类的有效数目

的计算公式为：

其中，N表示有效方位单元数，

表示类别矩阵L^(t-1)第i行的第j个元素，∑表示求和操作，i＝1，…，2N，j＝1，…，J^(t-1)，J^(t-1)表示干扰向量ε的聚类数；

混合系数π^(t)的第j个元素

的计算公式为：

其中，γ_j表示混合系数π的超参数，N_a表示缺损回波矩阵S_r的方位单元数；

S_r第q个距离单元散射点分布的权向量

的计算公式为：

其中，∑表示协方差矩阵，

表示第t-1次迭代的协方差期望，ψ＝[ψ₁，...，ψ_i，...，ψ_2N]表示过渡变量，且

表示精度向量α^(t-1)的第j个元素，

表示实傅里叶字典Φ的第i行，I表示单位矩阵，

表示实转置回波矩阵S的第f个列向量的第i个元素，

表示均值向量μ^(t-1)的第j个元素，[·]^-1表示矩阵求逆操作，T表示转置操作，f＝q；

均值向量μ^(t)的第j个元素

的计算公式为：

精度向量α^(t)的第j个元素

的计算公式为：

其中，a和b表示精度向量α第j个元素α_j的超参数；

类别矩阵L^(t)的第i行的第j个元素

的计算公式为：

其中，N(·)表示高斯分布的概率密度；

协方差期望

的第d个对角线元素

的计算公式为：

其中，

表示第q个距离单元散射点分布权向量

的第d个元素，P表示协方差矩阵∑第d个对角线元素η_d的超参数，且P的值为1/2N_a，κ_(·)(·)表示第二类修正贝塞尔函数，d＝1，...，2N_a；

(8d)更新类别矩阵L^(t)、干扰向量ε的聚类数J^(t)、混合系数π^(t)、均值向量μ^(t)和精度向量α^(t)，实现步骤为：

(8d1)以类别矩阵L^(t)第1，…，2N行中最大元素对应列的序号M₁，…，M_2N为元素构造标签向量[M₁，…，M_2N]^T；

(8d2)对比序号1，…，J^(t-1)与标签向量[M₁，…，M_2N]^T，得到r个序号l₁，…，l_r，其中，l₁，…，l_r表示序号1，…，J^(t-1)在标签向量[M₁，…，M_2N]^T中未出现过的序号；

(8d3)剔除类别矩阵L^(t)的第l₁，…，l_r个列向量以及混合系数π^(t)、均值向量μ^(t)和精度向量α^(t)的第l₁，…，l_r个元素，得到剔除后的类别矩阵L^(t)′、剔除后的混合系数π^(t)′、剔除后的均值向量μ^(t)′和剔除后的精度向量α^(t)′，将剔除后的类别矩阵L^(t)′的每个元素除以所在行元素之和，得到归一化剔除后的类别矩阵L^(t)″；

(8d4)更新类别矩阵L^(t)＝L^(t)″，更新干扰向量ε的聚类数J^(t)＝J^(t-1)-r，更新混合系数π^(t)＝π^(t)′，更新均值向量μ^(t)＝μ^(t)′，更新精度向量α^(t)＝α^(t)′；

(8e)判断当前迭代次数t＞1是否成立，若是，执行步骤(8f)；否则，令t＝t+1，并执行步骤(8c)：

(8f)判断

是否成立，若是，则S_r第q个距离单元散射点分布的权向量

并执行步骤(8g)；否则，令t＝t+1，并执行步骤(8c)，其中，||·||₂表示2-范数操作；

(8g)判断S_r当前距离单元序号q＝N_r是否成立，若是，得到S_r每个距离单元散射点分布的权向量

否则，令q＝q+1，并执行步骤(8b)；

步骤9)构造复权值矩阵W_c：

将S_r每个距离单元散射点分布的权向量

按行拼接，得到维数为2N_a×N_r的实权值矩阵W，并通过W构造复权值矩阵W_c：

W_c＝W₁+jW₂

其中，W₁表示实权值矩阵W第1，…，N_a行元素构成的矩阵，W₂表示实权值矩阵W第N_a+1，…，2N_a行元素构成的矩阵；

步骤10)获取运动目标的二维逆合成孔径雷达ISAR成像结果：

对复权值矩阵W_c进行转置，得到干扰环境下运动目标的二维高分辨逆合成孔径雷达ISAR成像结果。

以下结合仿真实验，对本发明的技术效果作进一步说明。

1.仿真条件和内容：

本发明的仿真实验采用工作在C波段的雷达，对应载频为10GHZ，带宽为0.4GHZ。目标长度为8.75米，翼展宽度为6米，回波数据的缺损率为50％。

仿真1，对距离脉冲压缩后的回波数据的第1～80列、207～320列以及451～512列产生缺损，并添加射频干扰，绘制距离脉冲压缩后的缺损回波的仿真图，其结果如图2(a)所示；

仿真2，利用现有技术中的正交匹配追踪OMP方法对缺损回波进行重构，绘制成像结果仿真图，其结果如图2(b)所示；

仿真3，利用现有技术中的相关向量机RVM对缺损回波进行重构，绘制成像结果仿真图，其结果如图2(c)所示；

仿真4，利用本发明对缺损回波进行重构，绘制成像结果仿真图，其结果如图2(d)所示。

2.仿真结果分析：

图2(a)为缺损率为50％的距离脉冲压缩后的缺损回波的仿真图，图2(a)中的横坐标表示距离脉冲压缩后的缺损回波的慢时间，纵坐标表示距离脉冲压缩后的缺损回波的距离单元。

图2(b)为利用现有技术中的正交匹配追踪OMP方法对缺损回波重构后的成像结果仿真图，图2(b)中的横坐标表示成像结果的方位单元，纵坐标表示成像结果的距离单元，可以看出利用现有的正交匹配追踪OMP方法得到的成像结果仿真图在第240到280列附近具有强干扰，虚假点较多，无法清楚地呈现目标的几何结构。

图2(c)为利用现有技术中的相关向量机RVM对缺损回波重构后的成像结果仿真图，图2(c)中的横坐标表示成像结果的方位单元，纵坐标表示成像结果的距离单元，可以看出利用现有的相关向量机RVM得到的成像结果仿真图受到残留干扰信号的影响，同样无法清楚地呈现目标的几何结构。

图2(d)为利用本发明对缺损回波重构后的成像结果仿真图，图2(d)中的横坐标表示成像结果的方位单元，纵坐标表示成像结果的距离单元，与图2(b)和图2(c)对比可得，相比于现有技术中的正交匹配追踪OMP方法和相关向量机RVM重构后的成像结果仿真图，利用本发明得到的重构后的成像结果仿真图能够清楚地呈现目标的几何结构，虚假点少，图像聚焦良好，分辨率显著提高。

由上述仿真结果表明，采用本发明得到的逆合成孔径雷达ISAR成像结果优于其它现有技术的逆合成孔径雷达ISAR成像结果。因此，与现有技术相比，本发明能够实现准确地环境统计建模及散射点先验的灵活建模，在强干扰环境下获得目标聚焦良好的ISAR像。

Claims (4)

1.一种基于环境统计建模的ISAR高分辨成像方法，其特征在于，包括如下步骤：

(1)获取运动目标的缺损回波矩阵S_r：

逆合成孔径雷达ISAR向运动目标发射线性调频信号，并接收干扰环境下运动目标维数为N_r×N_a的缺损回波矩阵S_r，其中，N_r表示缺损回波矩阵S_r的距离单元数，N_a表示缺损回波矩阵S_r的方位单元数；

(2)获取有效回波矩阵S_e：

沿距离向对缺损回波矩阵S_r进行脉冲压缩，得到距离脉冲压缩后的缺损回波矩阵S_pc，并剔除S_pc中缺损的列向量，得到维数为N_r×N的有效回波矩阵S_e，其中，N表示有效方位单元数；

(3)获取实转置回波矩阵S：

对有效回波矩阵S_e进行转置，得到复转置回波矩阵S_c，并通过S_c构造维数为2N×N_r的实转置回波矩阵S：

其中，Re(·)表示取实部操作，Im(·)表示取虚部操作；

(4)构造实傅里叶字典Φ：

(4a)以

为元素，构造维数为N_a×N_a的复傅里叶字典Φ_c，并剔除Φ_c中与缺损回波矩阵S_r的缺损列序号相对应的行，得到维数为N×N_a的有效傅里叶字典Φ_e，其中，e^(·)表示以自然常数为底的指数操作，j表示虚数单位，π表示圆周率，m表示复傅里叶字典Φ_c行的序号，n表示复傅里叶字典Φ_c列的序号，行序号m与列序号n的取值范围均为[-N_a/2,N_a/2-1]；

(4b)通过Φ_e构造实傅里叶字典Φ：

(5)构建缺损回波矩阵S_r每个距离单元的稀疏信号表征模型：

通过实傅里叶字典Φ和实转置回波矩阵S，构建缺损回波矩阵S_r中第q个距离单元的稀疏信号表征模型：

S_f＝Φω_q+ε

其中，S_f表示实转置回波矩阵S的第f个列向量，ω_q表示缺损回波矩阵S_r第q个距离单元散射点分布的权向量，ε表示干扰向量，f＝q＝1,…,N_r；

(6)构建干扰向量ε的高斯混合模型：

构建干扰向量ε的高斯混合模型，该模型包含干扰向量ε的第i个元素ε_i的概率分布p(ε_i|π,μ,α)、精度向量α的第j个元素α_j的先验分布p(α_j)和混合系数π的先验分布p(π)，其中：

干扰向量ε的第i个元素ε_i的概率分布p(ε_i|π,μ,α)的表达式为：

其中，J表示干扰向量ε的聚类数，π_j表示混合系数π的第j个元素，μ_j表示均值向量μ的第j个元素，α_j表示精度向量α的第j个元素，N(·)表示高斯分布的概率密度，Σ表示求和操作，i＝1,…,2N，j＝1,…,J；

精度向量α的第j个元素α_j的先验分布p(α_j)的表达式为：

p(α_j)＝Gam(α_j|a,b)

其中，a和b表示精度向量α的第j个元素α_j的超参数，Gam(·)表示伽马分布的概率密度；

混合系数π的先验分布p(π)的表达式为：

p(π)＝Dir(π|γ₁,...,γ_J)

其中，γ₁,...,γ_J表示混合系数π的超参数，Dir(·)表示狄利克雷分布的概率密度；

(7)构建权向量ω_q的伽马过程-高斯层级先验：

构建权向量ω_q的伽马过程-高斯层级先验，该先验包含权向量ω_q的先验分布p(ω_q)和协方差矩阵Σ第d个对角线元素η_d的先验分布p(η_d)，其中：

权向量ω_q的先验分布p(ω_q)的表达式为：

p(ω_q)＝N(ω_q|0,Σ)

其中，协方差矩阵Σ为对角阵；

协方差矩阵Σ第d个对角线元素η_d的先验分布p(η_d)的表达式为：

p(η_d)＝Gam(η_d|P,1)

其中，P表示协方差矩阵Σ第d个对角线元素η_d的超参数，d＝1,…,2N_a；

(8)利用最大后验-期望最大MAP-EM算法计算S_r每个距离单元散射点分布的权向量：

(8a)令S_r初始距离单元序号为q＝1；

(8b)令初始迭代次数为t＝1，并初始化参数：

在S_r第q个距离单元中，令干扰向量ε的聚类数J的初始值为10，令均值向量μ包含元素的个数为10，且各元素的初始值为随机值，令精度向量α包含元素的个数为10，且各元素的初始值均为0.01，令协方差矩阵Σ的初始值为对角线元素均为100的维数为2N_a×2N_a的对角阵，令混合系数π的超参数γ₁,...,γ_J均为1，令精度向量α每个元素的超参数a和b均为10^-4，令协方差矩阵Σ每个对角线元素的超参数P均为1/2N_a，令类别矩阵L的初始值的维数为2N×10，且类别矩阵L每行中仅有一个元素为1，其余元素均为0，每行非零元素的列序号在整数[1,10]中随机产生；

(8c)计算干扰向量ε第j个聚类的有效数目

混合系数π^(t)、S_r第q个距离单元散射点分布的权向量

均值向量μ^(t)、精度向量α^(t)、类别矩阵L^(t)和协方差期望

的值；

(8d)更新类别矩阵L^(t)、干扰向量ε的聚类数J^(t)、混合系数π^(t)、均值向量μ^(t)和精度向量α^(t)；

(8e)判断当前迭代次数t＞1是否成立，若是，执行步骤(8f)；否则，令t＝t+1，并执行步骤(8c)；

(8f)判断

是否成立，若是，则S_r第q个距离单元散射点分布的权向量

并执行步骤(8g)；否则，令t＝t+1，并执行步骤(8c)，其中，||·||₂表示2-范数操作；

(8g)判断S_r当前距离单元序号q＝N_r是否成立，若是，得到S_r每个距离单元散射点分布的权向量

否则，令q＝q+1，并执行步骤(8b)；

(9)构造复权值矩阵W_c：

将S_r每个距离单元散射点分布的权向量

按行拼接，得到维数为2N_a×N_r的实权值矩阵W，并通过W构造复权值矩阵W_c：

W_c＝W₁+jW₂

其中，W₁表示实权值矩阵W第1,…,N_a行元素构成的矩阵，W₂表示实权值矩阵W第N_a+1,…,2N_a行元素构成的矩阵；

(10)获取运动目标的二维逆合成孔径雷达ISAR成像结果：

对复权值矩阵W_c进行转置，得到干扰环境下运动目标的二维高分辨逆合成孔径雷达ISAR成像结果。

2.根据权利要求1所述基于环境统计建模的ISAR高分辨成像方法，其特征在于，步骤(2)中所述的距离脉冲压缩后的缺损回波矩阵S_pc，实现步骤为：

(2a)将逆合成孔径雷达ISAR至场景中心的距离作为参考距离，并选取载频和调频率与逆合成孔径雷达ISAR发射信号相同，且距离为参考距离的线性调频信号作为参考信号S_ref；

(2b)将参考信号S_ref取共轭后与接收的缺损回波矩阵S_r相乘，得到解线频调后的信号S_de，并对S_de沿快时间维作傅里叶变换，得到距离脉冲压缩后的缺损回波矩阵S_pc。

3.根据权利要求1所述基于环境统计建模的ISAR高分辨成像方法，其特征在于，步骤(8c)中所述的计算干扰向量ε第j个聚类的有效数目

混合系数π^(t)、S_r第q个距离单元散射点分布的权向量

均值向量μ^(t)、精度向量α^(t)、类别矩阵L^(t)和协方差期望

的值，计算公式分别为：

干扰向量ε第j个聚类的有效数目

的计算公式为：

其中,N表示有效方位单元数，

表示类别矩阵L^(t-1)第i行的第j个元素，Σ表示求和操作，i＝1,…,2N,j＝1,…,J^(t-1)，J^(t-1)表示干扰向量ε的聚类数；

混合系数π^(t)的第j个元素

的计算公式为：

其中，γ_j表示混合系数π的超参数，N_a表示缺损回波矩阵S_r的方位单元数；

S_r第q个距离单元散射点分布的权向量

的计算公式为：

其中，Σ表示协方差矩阵，

表示第t-1次迭代的协方差期望，ψ＝[ψ₁,...,ψ_i,...,ψ_2N]表示过渡变量，且

表示精度向量α^(t-1)的第j个元素，

表示实傅里叶字典Φ的第i行，I表示单位矩阵，

表示实转置回波矩阵S的第f个列向量的第i个元素，

表示均值向量μ^(t-1)的第j个元素，[·]^-1表示矩阵求逆操作，T表示转置操作，f＝q；

均值向量μ^(t)的第j个元素

的计算公式为：

精度向量α^(t)的第j个元素

的计算公式为：

其中，a和b表示精度向量α第j个元素α_j的超参数；

类别矩阵L^(t)的第i行的第j个元素

的计算公式为：

其中，N(·)表示高斯分布的概率密度；

协方差期望

的第d个对角线元素

的计算公式为：

其中，

表示第q个距离单元散射点分布权向量

的第d个元素，P表示协方差矩阵Σ第d个对角线元素η_d的超参数，且P的值为1/2N_a，κ_(·)(·)表示第二类修正贝塞尔函数，d＝1,...,2N_a。

4.根据权利要求1所述基于环境统计建模的ISAR高分辨成像方法，其特征在于，步骤(8d)中所述的更新类别矩阵L^(t)、干扰向量ε的聚类数J^(t)、混合系数π^(t)、均值向量μ^(t)和精度向量α^(t)，实现步骤为：

(8d1)以类别矩阵L^(t)第1,…,2N行中最大元素对应列的序号M₁,…,M_2N为元素构造标签向量[M₁,…,M_2N]^T；

(8d2)对比序号1,,J^(t-1)与标签向量[M₁,…,M_2N]^T，得到r个序号l₁,…,l_r，其中，l₁,…,l_r表示序号1,…,J^(t-1)在标签向量[M₁,…,M_2N]^T中未出现过的序号；

(8d3)剔除类别矩阵L^(t)的第l₁,…,l_r个列向量以及混合系数π^(t)、均值向量μ^(t)和精度向量α^(t)的第l₁,…,l_r个元素，得到剔除后的类别矩阵L^(t)′、剔除后的混合系数π^(t)′、剔除后的均值向量μ^(t)′和剔除后的精度向量α^(t)′，将剔除后的类别矩阵L^(t)′的每个元素除以所在行元素之和，得到归一化剔除后的类别矩阵L^(t)″；

(8d4)更新类别矩阵L^(t)＝L^(t)″，更新干扰向量ε的聚类数J^(t)＝J^(t-1)-r，更新混合系数π^(t)＝π^(t)′，更新均值向量μ^(t)＝μ^(t)′，更新精度向量α^(t)＝α^(t)′。

Images