CN114488025A - 基于平滑吕变换的空间目标三维成像方法 - Google Patents

Abstract

本发明提出了一种基于平滑吕变换的空间目标三维成像方法，用于解决现有技术中存在的成像结果依赖散射点关联的问题以及需要多个天线构成干涉式雷达的问题，实现步骤为：获取空间目标回波的时频矩向量；构建成像信号矩阵；获取成像信号矩阵的插值后RD像；对成像信号矩阵进行平滑吕变换；获取空间目标三维成像结果。本发明无需多个天线，且不需要关联散射点，得到空间目标准确的三维像结果。

Description

基于平滑吕变换的空间目标三维成像方法

技术领域

本发明属于雷达信号处理技术领域，特别涉及一种基于平滑吕变换的空间目标三维成像方法，可用于空间目标识别或分类。

背景技术

空间目标是指在大气层外围绕地球飞行的目标，常见的空间目标包含卫星，弹道目标，运载器碎片等。空间目标的识别技术是空间目标监视系统的重要环节，对维护空间安全有着重要的意义。对目标发射雷达信号，得到目标的雷达回波，使用三维成像算法处理雷达回波可以得到目标的三维图像矩阵，三维图像矩阵中的峰值点表征目标上散射点的位置，进一步可以估计目标尺寸特征，用于空间目标识别与分类。目前三维成像算法分为单雷达成像算法与多雷达成像算法，其中：单雷达成像算法构建三维成像空间，使用信号变换方法处理多个脉冲的雷达回波信号，得到三维成像矩阵，但是，该方类法受到空间目标运动影响，得到的成像结果模糊，难以提取出三维成像矩阵中的峰值，得到目标散射点的位置；多雷达成像算法使用多个雷达的单脉冲回波信号，根据信号的相位差得到目标散射点位置，合成目标三维成像结果，但是，该类方法需要多个雷达设备，对硬件的要求高。

Di Mo，Ning Wang和Guangzuo Li等人在2019年IEEE Geoscience and RemoteSensing Letters上发表文章“3-D Inverse Synthetic Aperture Ladar Imaging andScaling of Space Debris Based on the Fractional Fourier Transform”中提出一种基于分数阶傅里叶变换的空间碎片目标三维成像方法。该方法是一种单雷达成像算法，使用散射点距离、多普勒频率和调频率构建三维空间，使用分数阶傅里叶变换处理雷达回波信号，得到三维图像矩阵。该方法具体步骤是：第一步，建立空间碎片三维几何模型，得到空间碎片的回波信号模型；第二步，根据回波信号模型，设置合适的分数阶参数；第三步，使用分数阶傅里叶变换处理空间碎片雷达回波信号，得到三维图像矩阵。该方法计算复杂度低，且只使用单个雷达接收目标回波，降低硬件复杂度，但是，该方法仍然存在的不足之处是，分数阶傅里叶变换无法解决目标运动带来的成像模糊问题，导致三维图像矩阵模糊，影响后续散射点提取以及目标识别。

中国人民解放军空军工程大学在其申请的专利文献“一种空间高速运动目标干涉式三维成像方法”(专利申请号：CN201710212479.9，申请公布号：CN107024684A)中提出了一种空间快速运动目标干涉式三维成像方法。干涉式三维成像方法是一种多雷达成像算法，该方法使用三个雷达接收目标雷达回波信号，并使用三个回波的相位差估计目标散射点三维空间位置，合成目标三维成像结果。该方法具体步骤是：第一步，对三天线回波信号经过平动补偿，进行解线调频处理，构建含有未知速度的联合参数化稀疏表征模型；第二步，采用一种改进的OMP算法对联合参数化稀疏表征模型进行求解，获得运动速度、各散射点的距离向未知以及干涉相位信息；第三步，利用获得的干涉相位信息进行干涉处理，得到各散射点方位向和俯仰向位置，实现空间目标的三维成像。该方法只需要一个脉冲的雷达回波信号即可完成目标的三维成像，不受目标运动影响，但是，该方法仍然存在的不足之处是，需要三个天线的雷达接收目标回波信号，硬件复杂度高。

发明内容

本发明的目的是针对上述现有技术存在的不足，提出了一种基于平滑吕变换的空间目标三维成像方法，用于解决现有技术中存在的成像模糊问题以及需要多个天线构成干涉式雷达的技术问题。

实现本发明目的的思路是，选择合适的成像时间，分别使用RD成像算法和平滑吕变换方法成像，再使用Khatri-Rao积对成像结果作三维拼接处理，得到目标三维成像结果。其中，选择合适的成像时间，分别使用RD成像算法和平滑吕变换方法成像，再使用Khatri-Rao积对成像结果作三维拼接处理，具体步骤为：

(1)获取空间目标回波的时频矩向量S^Me：

(1a)对雷达发射的M个脉冲信号的回波信号各进行N次采样，得到维数为N×M的回波信号矩阵S，然后以w_STFT为短时傅里叶变换的窗长，以K为时频谱频率分辨单元数对S中的第N行各个元素进行短时傅里叶变换，得到维度为K×I的时频谱矩阵，并将每一列的时频谱组成时频谱向量

其中：N≥2，M≥2，1≤w_STFT≤M，K≥2，I为时频谱脉冲数，I＝M-w_STFT+1，

表示第i列第k个时频谱，i＝1,2,...,I，k＝1,2,...,K；

(1b)计算每个时频谱向量

的时频矩

得到长度为I的时频矩向量

(2)构建成像信号矩阵：

获取时频矩向量的最小时频矩索引α＝argmin(S^Me)，并在雷达发射的M个脉冲信号的回波信号中以α为初始点，连续选取M′个脉冲信号的回波信号作为成像信号，得到维数为N×M′的成像信号矩阵S′，其中，argmin(·)表示获取向量最小元素索引操作，2≤M′≤M-α；

(3)获取成像信号矩阵S′的插值后RD像：

对成像信号矩阵S′前W列的每个元素作RD成像处理，得到维数为D×W的RD像P，并对RD像每个像素值作线性插值处理，得到维数为D×W′的插值后的RD像P′，其中：D表示RD像的距离向单元数，D＝N，W＝ceil(M′/2)，W′＝M′；

(4)对成像信号矩阵S′进行平滑吕变换：

(4a)构建维数为X×Y、第x行第y列的元素为U_xy的尺度变换信号矩阵U，以及第x行第y列的元素为V_xy的尺度变换信号矩阵V，其中，U_xy和V_xy的表达式分别为：

其中：X＝M′，x＝1,2,...,X，Y＝M′，y＝1,2,...,Y，Z表示信号最小频点数，Z≥100，z＝1,2,...,Z，S′_Nm表示成像信号矩阵S′的第N行第m列元素，T_s表示回波脉冲间隔时间，(·)^*表示取共轭操作；

(4b)对尺度变换信号矩阵U和V的各个元素作以f_th为窗长的平滑卷积操作，得到维数为X×Y′尺度变换时间-中心频率像P^ST，其中：第x行第y′列像素值表达式为

其中：Y′＝M′，y′＝1,2,...,Y′，f_th为整数，1≤f_th≤M′且f_th≤5，

表示时域尺度变换信号矩阵U的第x行第y₁列的元素，y₁＝-y′+y，

表示时域尺度变换信号矩阵的第x行第y₂列的元素V，y₂＝y′-y；

(4c)对尺度变换时间-中心频率像P^ST各行作快速傅里叶变换，得到维数为X×Y′的中心频率-调频率图P″；

(5)获取空间目标三维成像结果：

对插值后的RD像P′的每个元素和中心频率-调频率图P″的每个元素作Khatri-Rao积，得到维数为N×M′×M′的三维图像T。

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

1)本发明对插值后的RD像的每个元素和通过对成像信号矩阵进行平滑吕变换所获取的中心频率-调频率图的每个元素作Khatri-Rao积，以获取空间目标三维成像结果，由于平滑吕变换使用了加窗平滑卷积操作，可以抑制三维图像矩阵峰值旁瓣，克服了现有技术中存在的由目标运动造成的成像模糊问题，使得本发明得到三维成像矩阵峰值能量更加集中；

2)本发明通过对单个雷达发射的M个脉冲信号的回波信号各进行N次采样，获得回波信号矩阵，雷达形式简单，硬件复杂度低，克服了现有技术中需要三个天线构成干涉式雷达的问题，使得本发明可以在只有单雷达天线的情况下，仍然可以得到准确的三维成像结果。

附图说明

图1是本发明的实现流程图。

图2是本发明实施例中插值后的RD成像结果图。

图3是本发明实施例中中心频率-调频率图结果图。

图4是本发明仿真实验采用的目标几何模型图。

图5是本发明仿真实验获得的三维成像结果图。

图6是现有技术仿真实验获得的三维成像结果图。

具体实施方式

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

参照图1，本发明包括如下步骤：

步骤1)获取空间目标回波的时频矩向量S^Me：

步骤1a)对雷达发射的M个脉冲信号的回波信号各进行N次采样，得到维数为N×M的回波信号矩阵S，然后以w_STFT为短时傅里叶变换的窗长，以K为时频谱频率分辨单元数对S中的第N行各个元素进行短时傅里叶变换，得到维度为K×I的时频谱矩阵，并将每一列的时频谱组成时频谱向量

其中：N≥2，M≥2，1≤w_STFT≤M，K≥2，I为时频谱脉冲数，I＝M-w_STFT+1，

表示第i列第k个时频谱，i＝1,2,...,I，k＝1,2,...,K，在本实施例中，M＝2000，N＝200，w_STFT＝32，K＝1024，本发明使用的是单个雷达M个脉冲的回波信号，雷达形式简单，硬件复杂度低，克服了现有技术中需要三个天线构成干涉式雷达的问题，使得本发明可以在只有单雷达天线的情况下，仍然可以得到准确的三维成像结果；

其中，第i列第k个时频谱

其计算公式为：

其中：exp(·)表示以自然常数为底的指数操作，j表示虚数单位，S_Nm表示回波信号矩阵S的第N行第m列元素，m＝1,2,...,M；

步骤1b)计算每个时频谱向量

的时频矩

得到长度为I的时频矩向量

其中，计算每个时频谱向量

的时频矩

计算公式为：

步骤2)构建成像信号矩阵：

获取时频矩向量的最小时频矩索引α＝argmin(S^Me)，并在雷达发射的M个脉冲信号的回波信号中以α为初始点，连续选取M′个脉冲信号的回波信号作为成像信号，得到维数为N×M′的成像信号矩阵S′，其中，argmin(·)表示获取向量最小元素索引操作，2≤M′≤M-α，在本实施例中，M′＝400，本发明通过计算时频矩，使用最小时频矩确定成像信号，当目标回波时频矩最小时，其散射点的运动速度最小，距离变化率也是最小的，可以减轻目标运动带来的图像模糊问题；

步骤3)获取成像信号矩阵S′的插值后RD像：

对成像信号矩阵S′前W列的每个元素作RD成像处理，得到维数为D×W的RD像P，并对RD像每个像素值作线性插值处理，得到维数为D×W′的插值后的RD像P′，其中：D表示RD像的距离向单元数，D＝N，W＝ceil(M′/2)，W′＝M′，在本发明中，使用一半成像信号的列向量做RD成像处理，可以减少成像信号的驻留时间，在短驻留时间内，目标运动导致的距离变化很小，减轻图像模糊，同时使用插值的方法，让RD像的维度恢复到原RD像大小；

其中，对成像信号矩阵S′前W列的每个元素作RD成像处理，第d行第w列像素值的RD成像公式为：

其中：d＝1,2,...,D，w＝1,2,...,W，S′_nm′表示成像信号矩阵S′的第n行第m′列元素；

对RD像每个像素值作线性插值处理，其中插值后的RD像P′第d行第w′列像素值的计算公式为：

其中：w′＝1,2,...,W′，

表示RD像P中第d行第w₁列像素值，w₁＝w′/2，

表示RD像P中第d行第w₂列像素值，

表示RD像P中第d行第w₃列像素值，

在本实施例中，图2为插值后的RD像结果，图中横轴为距离轴，纵轴为多普勒频率轴，图中白色亮点表示该图像的峰值，表征了散射点在该图像上的位置；

步骤4)对成像信号矩阵S′进行平滑吕变换：

步骤4a)构建维数为X×Y、第x行第y列的元素为U_xy的尺度变换信号矩阵U，以及第x行第y列的元素为V_xy的尺度变换信号矩阵V，其中，U_xy和V_xy的表达式分别为：

其中：X＝M′，x＝1,2,...,X，Y＝M′，y＝1,2,...,Y，Z表示信号最小频点数，Z≥100，z＝1,2,...,Z，S′_Nm表示成像信号矩阵S′的第N行第m列元素，T_s表示回波脉冲间隔时间，(·)^*表示取共轭操作，在本实施例中，Z＝500，T_s＝0.004；

步骤4b)对尺度变换信号矩阵U和V的各个元素作以f_th为窗长的平滑卷积操作，得到维数为X×Y′尺度变换时间-中心频率像P^ST，其中：第x行第y′列像素值表达式为

其中：Y′＝M′，y′＝1,2,...,Y′，f_th为整数，1≤f_th≤M′且f_th≤5，

表示时域尺度变换信号矩阵U的第x行第y₁列的元素，y₁＝-y′+y，

表示时域尺度变换信号矩阵的第x行第y₂列的元素V，y₂＝y′-y，在本实施例中，f_th＝5，本发明使用加窗平滑操作，可以有效滤除峰值附近的旁瓣，消除成像模糊；

步骤4c)对尺度变换时间-中心频率像P^ST各行作快速傅里叶变换，得到维数为X×Y′的中心频率-调频率图P″，在本实施例中，图3为中心频率-调频率图结果，图中横轴为调频率轴，纵轴为多普勒频率轴，图中白色亮点表示该图像的峰值，表征了散射点在该图像上的位置；

步骤5)获取空间目标三维成像结果：

对插值后的RD像P′的每个元素和中心频率-调频率图P″的每个元素作Khatri-Rao积，得到维数为N×M′×M′的三维图像T；

其中，对插值后的RD像P′的每个元素和中心频率-调频率图P″的每个元素作Khatri-Rao积，其中第d行第w′列第x个像素值的计算公式为：

T_dw′x＝P′_dw′×P″_xw′，

其中：P″_xw′表示中心频率-调频率图P″的第x行第w′列像素值。

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

1、仿真实验条件：

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

仿真实验的软件平台为：Windows 10操作系统，FEKO 7.0和Matlab 2017a。

仿真实验所使用的雷达回波信号由FEKO 7.0生成，目标的几何模型如图4所示，目标为一个四棱锥，在O-XYZ平面中，其四个散射点坐标分别为(0.5m,0m,0m)，(0m,-0.5m,0m)，(-0.5m,0.5m,0m)，(0m,0m,1m)。雷达回波信号的参数为：中心频率10GHz，脉冲重复频率500Hz，采样率2.4GHz，带宽4GHz。

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

对本发明的成像过程进行仿真，其结果如图5所示。

对现有的基于分数阶傅里叶变换的成像过程进行仿真，其结果如图6所示。

为了表征峰值能量集中情况，我们使用峰值能量集中度δ，其计算公式如下：

其中，P_peak表示三维成像矩阵峰值的和，P_total表示三维成像矩阵所有元素的和。峰值能量集中度越大，则该成像方法得到的三维图像峰值能量更加集中，图像越不模糊。

参照图5，图中以距离、多普勒频率和调频率为坐标轴构造三维空间，距离维和多普勒维构成一个面，调频率维构成多个层，构造目标三维图像，图中展示了本发明三维图像矩阵中0.15倍峰值的等值面。

参照图6，图中以距离、多普勒频率和调频率为坐标轴构造三维空间，距离维和多普勒维构成一个面，调频率维构成多个层，构造目标三维图像，图中展示了现有的基于分数阶傅里叶变换的三维图像矩阵中0.15倍峰值的等值面。

比较图5和图6，本发明的峰值能量集中度为0.012，基于分数阶傅里叶变换的三维成像方法峰值能量集中度为0.0011，可以看出，本发明可以有效解决成像模糊问题，得到的三维图像峰值能量更加集中。

Claims (6)

1.一种基于平滑吕变换的空间目标三维成像方法，其特征在于，包括如下步骤：

(1)获取空间目标回波的时频矩向量S^Me：

(1a)对雷达发射的M个脉冲信号的回波信号各进行N次采样，得到维数为N×M的回波信号矩阵S，然后以w_STFT为短时傅里叶变换的窗长，以K为时频谱频率分辨单元数对S中的第N行各个元素进行短时傅里叶变换，得到维度为K×I的时频谱矩阵，并将每一列的时频谱组成时频谱向量

其中：N≥2，M≥2，1≤w_STFT≤M，K≥2，I为时频谱脉冲数，I＝M-w_STFT+1，

表示第i列第k个时频谱，i＝1,2,...,I，k＝1,2,...,K；

(1b)计算每个时频谱向量

的时频矩

得到长度为I的时频矩向量

(2)构建成像信号矩阵：

获取时频矩向量的最小时频矩索引α＝argmin(S^Me)，并在雷达发射的M个脉冲信号的回波信号中以α为初始点，连续选取M′个脉冲信号的回波信号作为成像信号，得到维数为N×M′的成像信号矩阵S′，其中，argmin(·)表示获取向量最小元素索引操作，2≤M′≤M-α；

(3)获取成像信号矩阵S′的插值后RD像：

对成像信号矩阵S′前W列的每个元素作RD成像处理，得到维数为D×W的RD像P，并对RD像每个像素值作线性插值处理，得到维数为D×W′的插值后的RD像P′，其中：D表示RD像的距离向单元数，D＝N，W＝ceil(M′/2)，W′＝M′；

(4)对成像信号矩阵S′进行平滑吕变换：

(4a)构建维数为X×Y、第x行第y列的元素为U_xy的尺度变换信号矩阵U，以及第x行第y列的元素为V_xy的尺度变换信号矩阵V，其中，U_xy和V_xy的表达式分别为：

其中：X＝M′，x＝1,2,...,X，Y＝M′，y＝1,2,...,Y，Z表示信号最小频点数，Z≥100，z＝1,2,...,Z，S′_Nm表示成像信号矩阵S′的第N行第m列元素，T_s表示回波脉冲间隔时间，(·)^*表示取共轭操作；

(4b)对尺度变换信号矩阵U和V的各个元素作以f_th为窗长的平滑卷积操作，得到维数为X×Y′尺度变换时间-中心频率像P^ST，其中：第x行第y′列像素值表达式为

其中：Y′＝M′，y′＝1,2,...,Y′，f_th为整数，1≤f_th≤M′且f_th≤5，

表示时域尺度变换信号矩阵U的第x行第y₁列的元素，y₁＝-y′+y，

表示时域尺度变换信号矩阵的第x行第y₂列的元素V，y₂＝y′-y；

(4c)对尺度变换时间-中心频率像P^ST各行作快速傅里叶变换，得到维数为X×Y′的中心频率-调频率图P″；

(5)获取空间目标三维成像结果：

对插值后的RD像P′的每个元素和中心频率-调频率图P″的每个元素作Khatri-Rao积，得到维数为N×M′×M′的三维图像T。

2.根据权利要求1所述的基于平滑吕变换的空间目标三维成像方法，其特征在于，步骤(1a)中所述的第i列第k个时频谱

其计算公式为：

其中：exp(·)表示以自然常数为底的指数操作，j表示虚数单位，S_Nm表示回波信号矩阵S的第N行第m列元素，m＝1,2,...,M。

3.根据权利要求1所述的基于平滑吕变换的空间目标三维成像方法，其特征在于，步骤(1b)中所述的计算每个时频谱向量

的时频矩

计算公式：

4.根据权利要求1所述的基于平滑吕变换的空间目标三维成像方法，其特征在于，步骤(3)中所述的对成像信号矩阵S′前W列的每个元素作RD成像处理，其中第d行第w列像素值的RD成像公式为：

其中：d＝1,2,...,D，w＝1,2,...,W，S′_nm′表示成像信号矩阵S′的第n行第m′列元素。

5.根据权利要求1所述的基于平滑吕变换的空间目标三维成像方法，其特征在于，步骤(3)中所述的对RD像每个像素值作线性插值处理，其中插值后的RD像P′第d行第w′列像素值的计算公式为：

其中：w′＝1,2,...,W′，

表示RD像P中第d行第w₁列像素值，w₁＝w′/2，

表示RD像P中第d行第w₂列像素值，

表示RD像P中第d行第w₃列像素值，

6.根据权利要求1所述的基于平滑吕变换的空间目标三维成像方法，其特征在于，步骤(5)中所述的对插值后的RD像P′的每个元素和中心频率-调频率图P″的每个元素作Khatri-Rao积，其中第d行第w′列第x个像素值的计算公式为：

T_dw′x＝P′_dw′×P″_xw′，

其中：P″_xw′表示中心频率-调频率图P″的第x行第w′列像素值。

Images