CN113093140B - 基于加权admm的微动目标时频图重构方法 - Google Patents

Abstract

本发明提出了一种基于加权交替方向乘子法ADMM的微动目标时频图重构方法，解决现有技术中时频字典过大导致时频图重构速度慢，以及软阈值参数为定值导致重构图聚焦性下降的问题。本发明的实现步骤为：(1)获取微动目标的雷达回波信号；(2)截取雷达回波子信号；(3)生成时频字典；(4)生成每个子信号的加权向量；(5)计算每个子信号的时频分布向量；(6)构建微动目标的重构时频图。本发明生成了规模较小且唯一的时频字典，利用加权交替方向乘子法ADMM进行计算，提高了时频图的重构速度，并提高了重构图的聚焦性。

Description

基于加权ADMM的微动目标时频图重构方法

技术领域

本发明属于雷达技术领域，更进一步涉及雷达信号处理技术领域的一种基于加权交替方向乘子法ADMM(Alternating Direction Method of Multipliers)的微动目标时频图重构方法。本发明可用于对微动目标的时频图进行重构，得到聚焦性良好的时频图并将其用于微动目标的成像与微动目标的特征提取。

背景技术

微动是指目标或目标上的某些结构存在的独立于主体运动的振动或旋转，例如直升机叶片的旋转、弹道目标的自旋、进动和章动，以及人体行走时四肢的摆动等。微动会在目标主体运动对应的主多普勒周围产生边带频率调制，即产生微多普勒。对雷达回波信号进行短时傅里叶变换(STFT)后可以生成时频图，目前，常采用时频图对微多普勒进行描述。对微动目标进行时频重构获得的聚焦性良好的时频图，并将其用于微动目标的成像与微动目标时频图的特征提取。近年来，微动目标的雷达成像已受到雷达成像领域的广泛关注。

Chaowei Tong等人在其发表的论文“Ridge-Aware Weighted Sparse Time-Frequency Representation”(IEEE TRANSACTIONS ON SIGNAL PROCESSING,VOL.69,2021)中提出了一种基于脊感知加权的稀疏时频重构方法。该方法的具体步骤为：1.获取需要进行数值优化的信号；2.将回波加窗后的不同时刻的短时傅里叶变换矩阵堆叠成一个大矩阵，形成时频字典；3.对信号进行非线性时频变换，产生加权矩阵；4.运用快速迭代阈值收缩算法FISTA(Fast Iterative Shrinkage Threshold Algorithm)求解信号的时频分布。该方法存在的不足之处是，时频字典由回波加窗后的不同时刻的短时傅里叶变换矩阵堆叠而成，当雷达回波的脉冲数较多时，时频字典过大，同时FISTA算法迭代次数多，收敛速度慢，导致算法成本过高，重构速度慢。

西安电子科技大学在其申请的专利文献“基于时频字典的微动目标缺损回波高分辨成像方法”(专利申请号：201610034208.4，申请公开号：CN 105759264 A)中公开了一种基于时频字典的微动目标缺损回波高分辨成像方法。该方法的具体步骤为：1.对逆合成孔径雷达接收到的微动目标缺损回波数据进行距离压缩处理；2.产生回波的短时傅里叶变换矩阵；3.运用修正的增广拉格朗日算法计算每个距离单元缺损回波信号的时间-瞬时多普勒分布；4.将所有距离单元的时间-瞬时多普勒分布堆叠为三维矩阵，并选取不同时刻的二维矩阵生成目标的距离-瞬时多普勒图像序列。该方法存在的不足之处是，在使用修正的增广拉格朗日算法时使用的软阈值参数为定值，无法适应信号特点，容易过滤掉信号的特征信息，导致重构图的聚焦性下降。

发明内容

本发明的目的在于针对上述现有技术存在的不足，提出一种基于加权交替方向乘子法ADMM的微动目标时频图重构方法，用于解决现有技术中时频字典过大导致运算成本过高，重构速度慢，以及算法使用的软阈值参数为定值，无法适应信号特点，容易过滤掉信号特征信息，导致重构图的聚焦性下降的问题。

为实现上述目的，本发明的思路是：将雷达回波信号截取为多个雷达回波子信号，根据雷达回波子信号生成逆短时傅里叶变换矩阵，将其作为所有子信号共同的的时频字典，利用时频字典计算雷达回波中每个子信号对应的时频分布向量，以重构出时频图，该时频字典的规模较小，在计算过程中降低了运算复杂度，以提高重构速度；根据雷达回波子信号生成加权向量，将其作为软阈值参数，放大信号最大的时频系数同时抑制最小的时频系数，以提取信号的特征信息，提高重构图的聚焦性。

本发明的实现的具体步骤如下：

(1)获取微动目标的雷达回波信号：

获取雷达向微动目标发射脉冲信号后接收到的该微动目标的N×1维的雷达回波信号s₀，N表示雷达向微动目标发射的脉冲信号的总数；

(2)截取雷达回波子信号：

(2a)在雷达回波信号s₀前后各补(L-1)/2个零元素得到向量s，其中，L表示时间滑窗的窗长；

(2b)用窗长为L、滑动步长为1的时间滑窗在向量s上滑动，将每次滑动窗口中的所有元素组成该窗口的雷达回波向量；

(2c)在每个雷达回波向量后补N-L个零元素，得到该雷达回波的子信号；

(3)生成时频字典：

(3a)生成M×M维的离散傅里叶变换矩阵，其中，M与N的取值相等；

(3b)产生一个窗长为Len的窗函数向量，在该窗函数向量后补M-Len个零元素，以补零之后的窗函数向量中的元素为对角矩阵中对角元素得到一个对角矩阵，其中，Len与L的取值相等；

(3c)将离散傅里叶变换矩阵与对角矩阵进行矩阵进行相乘操作，得到短时傅里叶变换矩阵；

(3d)对短时傅里叶变换矩阵进行伪逆操作，得到逆短时傅里叶变换矩阵；

(4)生成每个子信号的加权向量：

(4a)使用窗函数向量对子信号进行短时傅里叶变换操作，得到子信号的时频向量；

(4b)使用窗函数的导数向量对子信号进行短时傅里叶变换操作，得到子信号的时频导数向量；

(4c)对子信号的时频导数向量和子信号的时频向量进行逐元素相除操作，得到子信号的非线性时频向量；

(4d)对非线性时频向量进行取绝对值操作，得到加权向量；

(5)计算每个子信号的时频分布向量：

(5a)利用分裂向量更新公式，计算当前迭代时的分裂向量；

(5b)利用辅助向量更新公式，计算当前迭代时的辅助向量；

(5c)将当前迭代时的分裂向量与辅助向量相加，得到当前迭代更新后的时频分布向量；

(5d)判断

的值是否小于终止阈值10^-3，若是，则终止迭代，执行步骤(6)，否则，将当前迭代次数加1后执行步骤(5a)；其中，||·||₂表示2-范数操作，f_k表示第k次迭代更新后的时频分布向量，/>

表示第k次迭代更新前的时频分布向量，当k＝0时，/>

为短时傅里叶变换矩阵与第i个雷达回波子信号的矩阵乘积；

(6)构建微动目标的重构时频图：

将所有更新后的雷达回波子信号的时频分布向量组成微动目标雷达回波信号s₀的重构时频图。

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

第一，由于本发明生成一个逆短时傅里叶变换矩阵作为时频字典，克服了现有技术中时频字典由回波加窗后的不同时刻的短时傅里叶变换矩阵堆叠而成，当雷达回波脉冲数较多时，时频字典大，运算成本高，重构速度慢的问题，使得本发明提高了时频图像的重建速度。

第二，本发明通过生成每个子信号的加权向量，克服了现有技术中算法使用的软阈值参数为定值，无法适应信号特点，容易过滤掉信号特征信息，导致重构图的聚焦性下降的问题，使得本发明能重构出聚焦性良好的时频图。

附图说明

图1为本发明的流程图；

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

具体实施方式

下面结合附图对本发明做进一步描述。

参照图1，对本发明的具体实施步骤做进一步描述。

步骤1，获取微动目标的雷达回波信号。

获取雷达向微动目标发射脉冲信号后接收到的该微动目标的N×1维的雷达回波信号s₀，N表示雷达向微动目标发射的脉冲信号的总数。

步骤2，截取雷达回波子信号。

在雷达回波信号s₀前后各补(L-1)/2个零元素得到向量s，其中，L表示时间滑窗的窗长。

用窗长为L、滑动步长为1的时间滑窗在向量s上滑动，将每次滑动窗口中的所有元素组成该窗口的雷达回波向量。

在每个雷达回波向量后补N-L个零元素，得到该雷达回波的子信号。

步骤3，生成时频字典。

利用下述公式计算M×M维的离散傅里叶变换矩阵中的各元素，其中，M与N的取值相等：

F(m,n)＝e^{-j2π(m-1)(n-1)/M}

其中，F(m,n)表示离散傅里叶变换矩阵F中的第m行、第n列的元素值，m,n∈[1,M]，e^(·)表示以自然指数e为底的指数操作，j表示虚数单位符号，π表示圆周率。

产生一个窗长为Len的窗函数向量，在该窗函数向量后补M-Len个零元素，以补零之后的窗函数向量中的元素为对角矩阵中对角元素得到一个对角矩阵，其中，Len与L的取值相等。

将离散傅里叶变换矩阵与对角矩阵进行矩阵进行相乘操作，得到短时傅里叶变换矩阵。

对短时傅里叶变换矩阵进行伪逆操作，得到逆短时傅里叶变换矩阵。

步骤4，生成每个子信号的加权向量。

使用窗函数向量对子信号进行短时傅里叶变换操作，得到子信号的时频向量。

使用窗函数的导数向量对子信号进行短时傅里叶变换操作，得到子信号的时频导数向量。

对子信号的时频导数向量和子信号的时频向量进行逐元素相除操作，得到子信号的非线性时频向量。

对非线性时频向量进行取绝对值操作，得到加权向量。

步骤5，计算每个子信号的时频分布向量。

5.1)利用下述分裂向量更新公式，计算当前迭代时的分裂向量：

其中，u_k表示第k次迭代更新后的分裂向量，soft(·)表示软阈值收缩操作，

表示第k次迭代更新前的辅助向量，λ表示取值为3的正则化系数，μ表示取值为30的辅助系数，w_i表示第i个加权向量，当k＝0时，u_k为M×1维的全零向量，/>

为M×1维的全零向量。

5.2)利用下述辅助向量更新公式，计算当前迭代时的辅助向量：

d_k＝G^H(s_i-Gu_k)/μ

其中，d_k表示第k次迭代更新后的辅助向量，G表示逆短时傅里叶变换矩阵，H表示共轭转置操作，s_i表示第i个雷达回波子信号。

5.3)将当前迭代时的分裂向量与辅助向量相加，得到当前迭代时的时频分布向量。

5.4)判断

的值是否小于终止阈值10^-3，若是，则终止迭代，执行步骤6，否则，将当前迭代次数加1后执行本步骤的5.1)；其中，||·||₂表示2-范数操作，f_k表示第k次迭代更新后的时频分布向量，/>

表示第k次迭代更新前的时频分布向量，当k＝0时，/>

为短时傅里叶变换矩阵与第i个雷达回波子信号的矩阵乘积。

步骤6，构建微动目标的重构时频图。

将所有更新后的雷达回波子信号的时频分布向量组成微动目标雷达回波信号s₀的重构时频图。

下面结合仿真实验对本发明的效果做进一步的说明：

1、仿真实验条件：

本发明的仿真实验的硬件平台为：处理器为Intel(R)Core(TM)i7-7700 CPU，主频为3.60GHz，内存8GB。

本发明的仿真实验的软件平台为：Windows 7操作系统和MATLAB2019a。

仿真实验均采用由MATLAB2019a软件仿真的信号作为雷达回波信号，雷达工作在C波段，对应载频为10GHZ，脉冲重复频率为512Hz，观测时间为0.5s，接收有四个散射点的自旋微动目标的雷达回波数据，雷达回波信号长度为256。

2、仿真内容及其结果分析：

本发明的仿真内容是分别采用本发明与现有技术的方法对仿真实验条件中产生的自旋微动目标的回波数据进行时频图重构。

在仿真实验中，采用的现有技术是指，西安电子科技大学在其申请的专利文献“基于时频字典的微动目标缺损回波高分辨成像方法”(专利申请号：201610034208.4，申请公开号：CN 105759264 A)中公开的一种基于时频字典的微动目标缺损回波高分辨成像方法。

本发明的仿真实验有2个。

仿真实验1：

利用现有技术中的基于时频字典的微动目标缺损回波高分辨成像方法，对自旋微动目标的雷达回波进行时频图重构，重构时间为1.14s，绘制其时频图重构结果，结果如图2(a)中白色线条组成的图像所示。

仿真实验2：

利用本发明对自旋微动目标的雷达回波进行时频图重构，重构时间为0.95s，绘制其时频图重构结果，结果如图2(b)中白色线条组成的图像所示。

下面结合图2的仿真图对本发明的效果做进一步的描述。

图2(a)中白色线条组成的图像为利用现有技术中的基于时频字典的微动目标缺损回波高分辨成像方法，对自旋微动目标的雷达回波进行重构的时频图。图2(a)中的横坐标表示时间滑窗单元，纵坐标表示多普勒单元。由图2(a)可以看出，利用现有技术中的基于时频字典的微动目标缺损回波高分辨成像方法得到的时频图像线条较粗，说明聚焦性不好。

图2(b)中白色线条组成的图像为利用本发明所提出的方法，对自旋微动目标的雷达回波进行重构的时频图。图2(b)中的横坐标表示时间滑窗单元，纵坐标表示多普勒单元。由图2(b)与图2(a)对比可见，利用本发明得到的时频图线条较细且清晰，聚焦性更好。

由上述仿真结果表明，本发明的方法生成逆短时傅里叶变换矩阵作为时频字典，在计算过程中提高了时频图重构速度，生成每个子信号的加权向量作为软阈值参数，能得到聚焦性更好的时频图。

Claims (4)

1.一种基于加权ADMM的微动目标时频图重构方法，其特征在于，生成时频字典，生成每个子信号的加权向量，计算雷达回波中每个子信号对应的时频分布向量；该方法的具体步骤如下：

(1)获取微动目标的雷达回波信号：

获取雷达向微动目标发射脉冲信号后接收到的该微动目标的N×1维的雷达回波信号s₀，N表示雷达向微动目标发射的脉冲信号的总数；

(2)截取雷达回波子信号：

(2a)在雷达回波信号s₀前后各补(L-1)/2个零元素得到向量s，其中，L表示时间滑窗的窗长；

(2b)用窗长为L、滑动步长为1的时间滑窗在向量s上滑动，将每次滑动窗口中的所有元素组成该窗口的雷达回波向量；

(2c)在每个雷达回波向量后补N-L个零元素，得到该雷达回波的子信号；

(3)生成时频字典：

(3a)生成M×M维的离散傅里叶变换矩阵，其中，M与N的取值相等；

(3b)产生一个窗长为Len的窗函数向量，在该窗函数向量后补M-Len个零元素，以补零之后的窗函数向量中的元素为对角矩阵中对角元素得到一个对角矩阵，其中，Len与L的取值相等；

(3c)将离散傅里叶变换矩阵与对角矩阵进行矩阵相乘操作，得到短时傅里叶变换矩阵；

(3d)对短时傅里叶变换矩阵进行伪逆操作，得到逆短时傅里叶变换矩阵；

(4)生成每个子信号的加权向量：

(4a)使用补零之后的窗函数向量对子信号进行短时傅里叶变换操作，得到子信号的时频向量；

(4b)使用窗函数的导数向量对子信号进行短时傅里叶变换操作，得到子信号的时频导数向量；

(4c)对子信号的时频导数向量和子信号的时频向量进行逐元素相除操作，得到子信号的非线性时频向量；

(4d)对非线性时频向量进行取绝对值操作，得到加权向量；

(5)计算每个子信号的时频分布向量：

(5a)利用分裂向量更新公式，计算当前迭代时的分裂向量；

(5b)利用辅助向量更新公式，计算当前迭代时的辅助向量；

(5c)将当前迭代时的分裂向量与辅助向量相加，得到当前迭代更新后的时频分布向量；

(5d)判断

的值是否小于终止阈值10^-3，若是，则终止迭代，执行步骤(6)，否则，将当前迭代次数加1后执行步骤(4a)；其中，||·||₂表示2-范数操作，f_k表示第k次迭代更新后的时频分布向量，/>

表示第k次迭代更新前的时频分布向量，当k＝0时，/>

为短时傅里叶变换矩阵与第i个雷达回波子信号的矩阵乘积；

(6)构建微动目标的重构时频图：

将所有更新后的雷达回波子信号的时频分布向量组成微动目标雷达回波信号s₀的重构时频图。

2.根据权利要求1所述的基于加权ADMM的微动目标时频图重构方法，其特征在于，步骤(3a)中所述的离散傅里叶变换矩阵中的各元素是由下式计算得到的：

F(m,n)＝e^{-j2π(m-1)(n-1)/M}

其中，F(m,n)表示离散傅里叶变换矩阵F中的第m行、第n列的元素值，m,n∈[1,M]，e^(·)表示以自然指数e为底的指数操作，j表示虚数单位符号，π表示圆周率。

3.根据权利要求1所述的基于加权ADMM的微动目标时频图重构方法，其特征在于，步骤(5a)中所述的分裂向量更新公式如下：

其中，u_k表示第k次迭代更新后的分裂向量，soft(·)表示软阈值收缩操作，

表示第k次迭代更新前的辅助向量，λ表示取值为3的正则化系数，μ表示取值为30的辅助系数，w_l表示第l个加权向量，l与i的取值对应相等，当k＝0时，u_k表示M×1维的全零向量，/>

表示M×1维的全零向量。

4.根据权利要求3所述的基于加权ADMM的微动目标时频图重构方法，其特征在于，步骤(5b)中所述的辅助向量更新公式如下：

d_k＝G^H(s_i-Gu_k)/μ

其中，d_k表示第k次迭代更新后的辅助向量，G表示逆短时傅里叶变换矩阵，H表示共轭转置操作，s_i表示第i个雷达回波子信号。

Images