CN105738894A - 基于增广拉普拉斯算子的微动群目标高分辨成像方法 - Google Patents

Abstract

本发明公开了一种基于增广拉普拉斯算子的微动群目标高分辨成像方法，主要解决现有的微动群目标高分成像方法对复杂微动形式不具有普适性，不适合在复杂环境下高分辨成像的问题，其实现步骤包括：(1)获取初始航迹矩阵；(2)利用增广拉普拉斯算子法重构初始航迹矩阵；(3)基于分类矩阵对多目标进行分类；(4)基于矩阵奇异值分解实现微动群目标高分辨成像。本发明具有稳定度高，成像误差小，对自旋、进动和章动普适性好的优点，可用于在初始航迹矩阵存在数据缺失，噪声干扰和奇异值这些复杂情况下的群微动目标高分辨成像。

Description

基于增广拉普拉斯算子的微动群目标高分辨成像方法

技术领域

本发明属于信号处理技术领域，更进一步涉及微动群目标成像中的基于增广拉普拉斯算子的微动群目标高分辨成像方法。本发明可以实现对多个弹道目标、空间碎片等微动目标在相互遮挡，噪声干扰等多种复杂环境下的高分辨成像。

背景技术

当采用高分辨雷达对微动目标成像时，在同一雷达波束宽度内可能包含多个微动目标，通常将这些同一雷达波束中具有相似运动状态的多个目标定义为群目标。

西安电子科技大学在其申请的专利“基于参数化的匀加速运动刚体群目标成像方法”(申请公布号：CN102778680A，申请号：201210198806.7)中提出了一种基于参数化的匀加速运动刚体微动群目标成像方法，该方法的具体步骤是，首先，雷达录取回波，并做模糊数估计与补偿以及二次相位粗补偿，其次做子目标中心的判定，图像的分割，以及相位项联合估计，最后实现群目标图像的合成。该方法运动补偿精确，图像聚焦良好，运算量较少，运算效率高。但是，仍然存在的不足之处是，这种群目标成像方法仅适合于匀加速运动的刚体群目标，对自旋、进动、章动等微动形式不具有普适性。

朱江等人在其发表的文献“基于块稀疏的空间碎片群目标成像方法”(电子信息学报，Vol.37No.3Mar.2015)中提出了一种基于块稀疏的群目标高分辨ISAR成像方法。该方法的具体步骤是，首先，基于块稀疏压缩感知理论，通过利用空间碎片群目标特性，抽取出各个碎片的高分辨一维距离像数据，其次结合平动补偿和距离多普勒算法得到各碎片的逆合成孔径雷达ISAR像。在空间碎片的数据可分离时，该方法能够有效实现群目标高分辨ISAR成像。但是，仍然存在的不足之处是，由于没有考虑群目标之间相互遮挡导致的回波数据缺失等复杂成像条件，使得其不能够在复杂环境下实现微动群目标的高分辨成像。

发明内容

本发明的目的在于克服上述现有技术的不足，提出一种基于增广拉普拉斯算子的微动群目标高分辨成像方法，以在回波数据缺失，噪声干扰，奇异值存在这些复杂环境下实现微动群目标的高分辨成像。

本发明的基本思路是：利用增广拉普拉斯算子法，重构在数据缺失，噪声干扰，奇异值存在的复杂环境下获取的初始航迹矩阵，并获得分类矩阵；基于分类矩阵构造图拉普拉斯算子，并用k均值方法将图拉普拉斯算子的特征向量聚类；通过带约束的矩阵奇异值分解的方法实现微动群目标的高分辨成像。其实现步骤包括如下：

(1)获取含有数据缺失，噪声、奇异值的初始航迹矩阵D，并记录该初始航迹矩阵中缺失数据的位置Ω；

(2)利用增广拉普拉斯算子法实现初始航迹矩阵D的重构；

(2a)令稀疏参数λ＝0.7，步长ρ＝1.1，精度δ＝10^-6；

(2b)令迭代次数k＝1，令第1次迭代的罚参数μ₁＝0.1，令第1次迭代的拉普拉斯乘法矩阵Λ₁、第1次迭代的奇异值矩阵E₁、第1次迭代的填补矩阵Z₁均为行数和列数与初始航迹矩阵D相同的零矩阵；

(2c)通过增广拉普拉斯算法循环迭代，获得第k次迭代的拉普拉斯乘法矩阵Λ_k，第k次迭代的奇异值矩阵E_k，第k次迭代的填补矩阵Z_k，最终得到重构矩阵A和分类矩阵C；

(3)基于分类矩阵的多微动目标分类：

(3a)将分类矩阵C的每一个元素取绝对值得到仿射矩阵M；

(3b)构造图拉普拉斯算子L并由图拉普拉斯算子零特征值的个数确定微动目标的个数；

(3c)将图拉普拉斯算子的特征向量用k均值方法聚类，得到每一个散射点对应的类别；

(3d)将属于同一类别的散射点航迹合并获得每个微动目标对应的子航迹矩阵；

(4)基于子航迹矩阵奇异值分解实现微动群目标的高分辨成像：

(4a)令微动目标序号j＝1；

(4b)将第j个子航迹矩阵中每个元素都减去该元素所在行所有元素的均值，得到平动校正后的第j个子航迹矩阵；

(4c)对平动校正后的第j个子航迹矩阵进行带约束的矩阵奇异值分解，得到第j个微动目标的散射点分布矩阵和雷达视线矩阵；

(4d)将第j个微动目标的散射点分布矩阵作为第j个微动目标的高分辨图像；

(4e)判断是否得到所有微动目标的高分辨图像，若是，执行步骤(4f)，否则，j＝j+1，返回步骤(4b)；

(4f)将每个微动目标的高分辨图像合并得到微动群目标的高分辨图像。

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

第一，本发明通过带约束的矩阵分解实现微动目标高分辨成像，克服了参数化的匀加速运动刚体群目标成像方法仅适合于刚体匀加速运动这一单一运动形式的缺点，对自旋、进动、章动等复杂微动形式具有普适性；

第二，本发明通过增广拉普拉斯算子的方法实现初始航迹矩阵的重构，使得该方法适合在数据缺失，噪声干扰，奇异值存在这些复杂环境下的微动群目标高分辨成像，克服了块稀疏的空间碎片群目标成像方法中不能够在复杂环境下实现微动群目标高分辨成像的缺点。

附图说明

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

图2为本发明的仿真结果图。

具体实施方式

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

参照图1，本发明的实现步骤如下：

步骤1，获取含有数据缺失，噪声、奇异值的初始航迹矩阵D，并记录初始航迹矩阵中缺失数据的位置Ω。

步骤2，利用增广拉普拉斯算子法实现初始航迹矩阵的重构：

(2a)令稀疏参数λ＝0.7，步长ρ＝1.1，精度δ＝10^-6；

(2b)令迭代次数k＝1，令第1次迭代的罚参数μ₁＝0.1，令第1次迭代的拉普拉斯乘法矩阵Λ₁、第1次迭代的奇异值矩阵E₁、第1次迭代的填补矩阵Z₁均为行数和列数与初始航迹矩阵D相同的零矩阵；

(2c)通过增广拉普拉斯算子法循环迭代，获得第k次迭代的拉普拉斯乘法矩阵Λ_k，第k次迭代的奇异值矩阵E_k，第k次迭代的填补矩阵Z_k，最终得到重构矩阵A和分类矩阵C：

(2c1)按照下式求第k次迭代的酉矩阵U_k，第k次迭代的特征值矩阵Σ_k，第k次迭代的特征向量矩阵V_k，

$\{U_k,{\mathrm{\Σ}}_k,V_k\}=SVD(D_k+\frac{1}{{\mathrm{\μ}}_k}{\mathrm{\Λ}}_k-E_k-Z_k),$

其中，SVD(·)表示矩阵的特征值分解操作，D表示初始航迹矩阵，μ_k表示第k次迭代的罚参数，Λ_k表示第k次迭代的拉普拉斯乘法矩阵，E_k表示第k迭次代的奇异值矩阵，Z_k表示第k次迭代的填补矩阵；

矩阵的分解操作方法通常有，矩阵奇异值分解，矩阵特征值分解，矩阵满秩分解等，本实例在步骤(2c1)中使用矩阵的特征值分解方法；

(2c2)取出特征向量矩阵V_k的前i列，构成第k次迭代的子特征向量矩阵V1_k，其中，i表示第k次迭代的特征值矩阵Σ_k中特征值大于的个数；

(2c3)根据子特征向量矩阵V1_k计算归一化向量n；

$n=\frac{(I-V{1}_{k}V{1}_k^T)\·ones}{\left|\right|(I-V{1}_{k}V{1}_k^T)\·ones|{|}_2}$

其中，I表示单位矩阵，ones表示元素个数为初始航迹矩阵列数，值全为1的向量，T表示矩阵转至运算，||·||₂表示求向量的2范数操作；

(2c4)按照下式计算第k次迭代的重构矩阵A_k：

$A_k=U_k{H}_{\sqrt{\frac{2}{{\mathrm{\μ}}_k}}}\[{\mathrm{\Σ}}_k\]{V_k}^T,$

其中，H[·]表示阈值算子操作，表示利用阈值算子，将第k次迭代的特征值矩阵Σ_k中特征值小于阈值的值赋为0时得到的阈值特征值矩阵；

(2c5)根据按照下式计算第k次迭代的分类矩阵C_k；

C_k＝[V1_k,n][V1_k,n]^T；

(2c6)按照下式计算第k次迭代的更新矩阵R_k；

$R_k=D+\frac{1}{{\mathrm{\μ}}_k}{\mathrm{\Λ}}_k-A_k-E_k;$

(2c7)判断更新矩阵R_k的F范数是否小于精度δ：若是，则将第k+1次迭代的填充矩阵Z_k+1赋值为更新矩阵R_k的值；否则，将第k+1次迭代的填充矩阵Z_k+1中与初始航迹矩阵D中非缺失数据位置相同处元素的值赋为更新矩阵R_k对应位置元素的值，将数据缺失位置Ω元素的值赋为其中，Z_k+1(i,j)表示第k+1次迭代的填充矩阵Z_k+1第i行，第j列的值，R_k(i,j)表示第k次迭代的更新矩阵R_k的第i行，第j列的值；

更新矩阵R_k的F范数按下式计算：

$\left|\right|R_k|{|}_F=\sqrt{\underset{i}{\Σ}\underset{j}{\Σ}{\left(r_{ij}\right)}^2},$

其中，||·||_F表示求矩阵F范数的操作，Σ表示求和操作，r_ij表示更新矩阵R_k的第i行第j列的值；

(2c8)按照下式计算第k次迭代的差值矩阵Y_k：

$Y_k=D+\frac{1}{{\mathrm{\μ}}_k}{\mathrm{\Λ}}_k-A_k-Z_{k+1};$

(2c9)依次判断差值矩阵Y_k的每一列构成向量的2范数是否小于若是，则将第k+1次迭代的奇异值矩阵E_k+1中的此列元素的值赋为0；否则，将E_k+1此列的每个值赋为其中，E_k+1(α,β)表示第k+1次迭代的奇异值矩阵E_k+1第β列的第α个值，Y_k(:,β)表示第k次迭代的差值矩阵Y_k的第β列，Y_k(α,β)表示第k次迭代的差值矩阵Y_k的第β列的第α个值；

(2c10)按照下式计算第k+1次迭代的拉普拉斯乘法矩阵：

Λ_k+1＝Λ_k+μ_k(D-A_k-E_k+1-Z_k+1)；

(2c11)将步长ρ与第k次迭代的罚参数μ_k相乘，得到第k+1次迭代的罚参数μ_k+1；

(2c12)判断迭代过程是否已经收敛，若是，则将第k次迭代的重构矩阵A_k作为重构矩阵A，将第k次迭代的分类矩阵C_k作为分类矩阵C，执行步骤(3)，否则，令k＝k+1，返回步骤(2c1)。

步骤3，基于分类矩阵的多微动目标分类。

(3a)将分类矩阵C的每一个元素取绝对值，得到仿射矩阵M；

(3b)构造图拉普拉斯算子L，并由图拉普拉斯算子零特征值的个数确定微动目标的个数，

图拉普拉斯算子L的构造通过下式进行：

L＝I-N^-1/2MN^1/2，

其中，I表示单位矩阵，M表示仿射矩阵；N表示对角矩阵，它的每个对角线元素等于仿射矩阵M对应行的所有元素之和；

(3c)对图拉普拉斯算子的特征向量进行聚类，得到每一个散射点对应的类别；

传统的聚类方法有c均值聚类方法和k均值聚类方法，在本实例中使用k均值聚类方法，其具体步骤如下：

(3c1)随机选择l个图拉普拉斯算子的特征向量作为聚类中心，其中l表示微动目标个数；

(3c2)依次计算所有特征向量与l个聚类中心的欧氏距离，并将其与之欧氏距离最近的一个聚类中心聚为一类；

(3c3)将(3c2)中得到的每一类特征向量的均值作为新的聚类中心；

(3c4)判断新的聚类中心是否发生改变，若是，则将每一特征向量对应的类别作为对应散射点所属的类别，否则，返回(3c2)；

(3d)将属于同一类别的散射点航迹合并获得每个微动目标对应的子航迹矩阵；

(3e)将属于同一类别的散射点航迹合并，获得每个微动目标对应的子航迹矩阵。

步骤4，基于子航迹矩阵奇异值分解实现微动群目标的高分辨成像。

(4a)令微动目标序号j＝1；

(4b)将第j个子航迹矩阵中每个元素都减去该元素所在行所有元素的均值，得到平动校正后的第j个子航迹矩阵；

(4c)对平动校正后的第j个子航迹矩阵进行带约束的矩阵奇异值分解，得到第j个微动目标的散射点分布矩阵和雷达视线矩阵；

(4d)将第j个微动目标的散射点分布矩阵作为第j个微动目标的高分辨图像；

(4e)判断是否得到所有微动目标的高分辨图像，若是，执行步骤(4f)，否则，令j＝j+1，返回步骤(4b)；

(4f)将每个微动目标的高分辨图像进行合并，得到微动群目标的高分辨图像。

本发明的效果可通过一下仿真进一步说明：

1、仿真实验条件：

仿真在MATLAB7.0软件下进行，仿真数据的参数如下：

雷达带宽为2GHz，载频为10GHz，PRF为100Hz，观测时间为1s，方位向回波为100次，雷达视线的方位角为50.6°，俯仰角为45.3°，仿真数据中包含两个微动目标，目标1包含9个强散射点，微动形式为进动，目标2包含5个强散射点，微动形式为自旋。

2、仿真实验内容：

采用本发明对两个微动目标进行高分辨成像，结果如图2，其中：

图2(a)为仿真过程中所假设的两个微动目标散的射点分布图，其中目标1散射点所在位置用黑点表示，目标2散射点所在位置用星号表示；

图2(b)为仿真过程中在数据缺失，噪声干扰，奇异值存在这些复杂环境下获取的图2(a)所假设的两个微动目标散射点的初始航迹矩阵图；

图2(c)为仿真过程中通过增广拉普拉斯算子循环迭代的方法，实现对图2(b)所获取的初始航迹矩阵的重构结果图；

图2(d)为仿真过程中利用增广拉普拉斯算子法，完成对目标1的高分辨成像结果图；

图2(e)为仿真过程中利用增广拉普拉斯算子法，完成对目标2的高分辨成像结果图；

从图2可以看出，采用增广的拉普拉斯算子法能够实现对微动群目标的高分辨成像。

Claims (4)

1.基于增广拉普拉斯算子的微动群目标高分辨成像方法，包括：

(1)获取含有数据缺失，噪声、奇异值的初始航迹矩阵D，并记录该初始航迹矩阵中缺失数据的位置Ω；

(2)利用增广拉普拉斯算子法实现初始航迹矩阵D的重构；

(2a)令稀疏参数λ＝0.7，步长ρ＝1.1，精度δ＝10^-6；

(2b)令迭代次数k＝1，令第1次迭代的罚参数μ₁＝0.1，令第1次迭代的拉普拉斯乘法矩阵Λ₁、第1次迭代的奇异值矩阵E₁、第1次迭代的填补矩阵Z₁均为行数和列数与初始航迹矩阵D相同的零矩阵；

(2c)通过增广拉普拉斯算子法循环迭代，获得第k次迭代的拉普拉斯乘法矩阵Λ_k，第k次迭代的奇异值矩阵E_k，第k次迭代的填补矩阵Z_k，最终得到重构矩阵A和分类矩阵C；

(3)基于分类矩阵的多微动目标分类：

(3a)将分类矩阵C的每一个元素取绝对值得到仿射矩阵M；

(3b)构造图拉普拉斯算子L并由图拉普拉斯算子零特征值的个数确定微动目标的个数；

(3c)将图拉普拉斯算子的特征向量用k均值方法聚类，得到每一个散射点对应的类别；

(3d)将属于同一类别的散射点航迹合并获得每个微动目标对应的子航迹矩阵；

(4)基于子航迹矩阵奇异值分解实现微动群目标的高分辨成像：

(4a)令微动目标序号j＝1；

(4b)将第j个子航迹矩阵中每个元素都减去该元素所在行所有元素的均值，得到平动校正后的第j个子航迹矩阵；

(4c)对平动校正后的第j个子航迹矩阵进行带约束的矩阵奇异值分解，得到第j个微动目标的散射点分布矩阵和雷达视线矩阵；

(4d)将第j个微动目标的散射点分布矩阵作为第j个微动目标的高分辨图像；

(4e)判断是否得到所有微动目标的高分辨图像，若是，执行步骤(4f)，否则，j＝j+1，返回步骤(4b)；

(4f)将每个微动目标的高分辨图像合并得到群微动目标的高分辨图像。

2.根据权利要求1所述的基于增广拉普拉斯算子的微动群目标高分辨成像方法

其特征在于，步骤(2c)中的循环迭代过程如下：

(2c1)按照下式求第k次迭代的酉矩阵U_k，第k次迭代的特征值矩阵Σ_k，第k次迭代的特征向量矩阵V_k，

$\{U_k,{\mathrm{\Σ}}_k,V_k\}=SVD(D_k+\frac{1}{{\mathrm{\μ}}_k}{\mathrm{\Λ}}_k-E_k-Z_k),$

其中，SVD(·)表示矩阵的特征值分解操作，D表示初始航迹矩阵，μ_k表示第k次迭代的罚参数，Λ_k表示第k次迭代的拉普拉斯乘法矩阵，E_k表示第k迭次代的奇异值矩阵，Z_k表示第k次迭代的填补矩阵；

(2c2)取出特征向量矩阵V_k的前i列，构成第k次迭代的子特征向量矩阵V1_k，其中，i表示第k次迭代的特征值矩阵Σ_k中特征值大于的个数；

(2c3)根据子特征向量矩阵V1_k计算归一化向量n；

$n=\frac{(I-V{1}_{k}V{1}_k^T)\·ones}{\left|\right|(I-V{1}_{k}V{1}_k^T)\·ones|{|}_2}$

其中，I表示单位矩阵，ones表示元素个数为初始航迹矩阵列数，其值全为1的向量，T表示矩阵转至运算，||·||₂表示求向量的2范数操作；

(2c4)按照下式计算第k次迭代的重构矩阵A_k：

$A_k=U_k{H}_{\sqrt{\frac{2}{{\mathrm{\μ}}_k}}}\[{\mathrm{\Σ}}_k\]{V_k}^T$

其中，H[·]表示阈值算子操作，表示利用阈值算子，将第k次迭代的特征值矩阵Σ_k中特征值小于阈值的值赋为0时得到的阈值特征值矩阵；

(2c5)根据按照下式计算第k次迭代的分类矩阵C_k；

C_k＝[V1_k,n][V1_k,n]^T；

(2c6)按照下式计算第k次迭代的更新矩阵R_k；

$R_k=D+\frac{1}{{\mathrm{\μ}}_k}{\mathrm{\Λ}}_k-A_k-E_k;$

(2c7)判断更新矩阵R_k的F范数是否小于精度δ：若是，则将第k+1次迭代的填充矩阵Z_k+1赋值为更新矩阵R_k的值；否则，将第k+1次迭代的填充矩阵Z_k+1中与初始航迹矩阵D中非缺失数据位置相同处元素的值赋为更新矩阵R_k对应位置元素的值，将数据缺失位置Ω元素的值赋为其中，Z_k+1(i,j)表示第k+1次迭代的填充矩阵Z_k+1第i行，第j列的值，R_k(i,j)表示第k次迭代的更新矩阵R_k的第i行，第j列的值；

(2c8)按照下式计算第k次迭代的差值矩阵Y_k；

$Y_k=D+\frac{1}{{\mathrm{\μ}}_k}{\mathrm{\Λ}}_k-A_k-Z_{k+1};$

(2c9)依次判断差值矩阵Y_k的每一列构成向量的2范数是否小于若是，则将第k+1迭次代的奇异值矩阵E_k+1中的此列元素的值赋为0；否则，将E_k+1此列的每个值赋为其中，E_k+1(α,β)表示第次k+1迭代的奇异值矩阵E_k+1第β列的第α个值，Y_k(:,β)表示第k次迭代的差值矩阵Y_k的第β列，Y_k(α,β)表示第k次迭代的差值矩阵Y_k的第β列的第α个值。

(2c10)按照下式计算第k+1次迭代的拉普拉斯乘法矩阵：

Λ_k+1＝Λ_k+μ_k(D-A_k-E_k+1-Z_k+1)；

(2c11)将步长ρ与第k次迭代的罚参数μ_k相乘，得到第k+1次迭代的罚参数μ_k+1；

(2c12)判断迭代过程是否已经收敛，若是，将第k次迭代的重构矩阵A_k作为重构矩阵A，将第k次迭代的分类矩阵C_k分类矩阵C，执行步骤(3)，否则，令k＝k+1，返回步骤(2c1)。

3.根据权利要求1所述的基于增广拉普拉斯算子的微动群目标高分辨成像方法

其特征在于，步骤(3b)中构造图拉普拉斯算子L，通过如下公式进行：

L＝I-N^-1/2MN^1/2

其中，I表示单位矩阵，M表示仿射矩阵；N表示对角矩阵，它的每个对角线元素等于仿射矩阵M对应行的所有元素之和。

4.根据权利要求1所述的基于增广拉普拉斯算子的微动群目标高分辨成像方法，其特征在于，步骤(3c)中用k均值方法对图拉普拉斯算子的特征向量进行聚类，按如下步骤进行：

(3c1)随机选择l个图拉普拉斯算子的特征向量作为聚类中心，其中l表示微动目标个数。

(3c2)依次计算所有特征向量与l个聚类中心的欧氏距离，并将其与之欧氏距离最近的一个聚类中心聚为一类；

(3c3)将(3c2)中得到的每一类特征向量的均值作为新的聚类中心；

(3c4)判断新的聚类中心是否发生改变，若是，则将每一特征向量对应的类别作为对应散射点所属的类别，否则，返回(3c2)。