CN114488058B - 一种基于stft谱图滑窗相消的微动杂波去除方法 - Google Patents

Abstract

本发明提出了一种基于STFT谱图滑窗相消的微动杂波去除方法。将时域数据进行短时傅里叶变换，得到STFT数据及STFT谱图；将STFT数据沿时间维等分为前、后参考窗数据，将STFT谱图沿时间维等分为前、后参考窗谱图；从STFT谱图截取多个滑窗谱图；将前、后参考窗谱图分别与各滑窗谱图相减，并计算相减后前、后参考窗谱图与原前、后参考窗谱图的强度比；通过设置门限确定前、后参考窗谱图中微动杂波位置；将微动杂波位置上的前、后参考窗数据置零以去除微动杂波；将去除微动杂波后的前、后参考窗数据分别沿时间维叠加，然后由逆傅里叶变换得到去除微动杂波后的前、后参考窗时域数据，将两者进行组合获取最终结果。本发明提高了雷达的探测性能及环境适应能力。

Description

一种基于STFT谱图滑窗相消的微动杂波去除方法

技术领域

本发明涉及雷达信号处理领域，特别是涉及一种基于STFT谱图滑窗相消的微动杂波去除方法。

背景技术

微多普勒效应(简称微动效应)是由目标或其部件的微动(如转动、振动、进动等)引起的频率调制现象。对雷达系统过来说，微动效应是一把“双刃剑”，一方面它能够反映目标的几何结构和运动状态，为目标参数估计、分类识别等应用提供宝贵的特征信息。另一方面，它也是雷达系统中的一种特殊杂波源，会抬高噪声基底，掩盖弱目标，造成虚警和漏检。

微动杂波是一种时变非平稳杂波，往往具有较大的多普勒拓展，常规动目标指示或时域杂波对消类方法很难取得令人满意的效果，而空域杂波抑制方法要求杂波和目标在角度维上是可分的，否则目标和杂波会被同时抑制。目前，微动杂波抑制方法主要分为参数化和非参数化方法。参数化方法要求对微动杂波进行精细建模，若模型无法反应杂波真实特性，将会因为模型失配导致算法性能下降甚至失效，因此参数化方法往往适用的条件较为苛刻，这在一定程度上限制了它的应用场景。相比较而言，非参数化方法对信号模型的依赖程度大大减弱，因此往往具有操作简便、适用范围广等优点。时频变换算法是一种常用的非平稳信号分析方法，因此被广泛应用于微动回波特性的分析中。微动杂波和目标回波在时频域中的形态表现具有明显差异。具体地，当积累时间较短时，可以认为目标回波在时频图中表现为在特定频率单元上的平行于时间轴的直线型能量条带，而微动杂波则会在一定频率范围内的多个频率单元上呈现复杂的分布形态。有鉴于此，在时频域中完成目标回波和微动杂波的分离，进而实现微动杂波去除的思路具有极大的发展前景。

发明内容

在详细分析目标回波和微动杂波在时频域中的形态差异，及现有时频域微动杂波去除方法的优缺点的基础上，本发明提出了一种基于STFT谱图滑窗相消的微动杂波去除方法，具体包含以下步骤：

步骤1：将原始时域采样数据进行短时傅里叶变换，得到STFT数据及STFT谱图；

步骤2：将STFT数据以时间维中心为对称轴等分为前参考窗数据、后参考窗数据，将STFT谱图以时间维中心为对称轴等分为前参考窗谱图、后参考窗谱图；

步骤3：设置滑窗数目K，从STFT谱图中截取得到K个滑窗谱图；

步骤4：将前参考窗谱图与各滑窗谱图进行对应元素相减得到相减后前参考窗谱图，计算相减后前参考窗谱图与原前参考窗谱图对应元素的比值；将后参考窗谱图与各滑窗谱图进行对应元素相减得到相减后后参考窗谱图，计算相减后后参考窗谱图与原后参考窗谱图对应元素的比值；结合设置的门限确定前参考窗谱图中微动杂波位置；结合设置的门限确定后参考窗谱图中微动杂波位置；

步骤5：将微动杂波所在位置上的前参考窗数据置零，得到去除微动杂波后的前参考窗数据；将微动杂波所在位置上的后参考窗数据置零，得到去除微动杂波后的后参考窗数据；

步骤6：将去除微动杂波后的前参考窗数据沿时间轴进行叠加得到前参考窗频谱，对前参考窗频谱进行逆傅里叶变换得到去除微动杂波后的前参考窗时域数据；将去除微动杂波后的后参考窗数据沿时间轴进行叠加得到后参考窗频谱，对后参考窗频谱进行逆傅里叶变换得到去除微动杂波后的后参考窗时域数据；将去除微动杂波后的前参考窗时域数据和去除微动杂波后的后参考窗时域数据进行组合得到去除微动杂波后的完整时域数据。

作为优选，步骤1所述将原始时域采样数据进行短时傅里叶变换为：

其中，N表示短时傅里叶变换的频率维数，M表示短时傅里叶变换的时间维数，S表示STFT数据，是一个N×M维的矩阵，S[m,k]表示第m行第k列STFT数据，x表示时域采样数据序列，所述时域采样数据矩阵的长度为N，x[n]表示第n个时域采样点，w表示窗函数序列，所述窗函数的有效长度为N_w，并通过补零将其拓展成长度为N的向量，它满足为常数的条件，w[n-mR]表示窗函数中的第n-mR个元素，R表示短时傅里叶变换相邻窗函数之间的步进长度；

将S表示为S＝[s₀,...,s_M-1]，s_m表示第m列SFFT数据；

通过Y[m,k]＝|S[m,k]|²,m∈[0,M-1],k∈[0,N-1]计算STFT谱图；

其中，Y表示STFT谱图，Y[m,k]表示STFT谱图在第m行第k列上的元素；

将Y表示为Y＝[y₀,...,y_M-1]，y_m表示第m列STFT谱图；

作为优选，步骤2所述前参考窗数据为：

S_F＝[s₀,...,s_M/2-1]

步骤2所述后参考窗数据为：

S_B＝[s_M/2,...,s_M-1]

步骤2所述前参考窗谱图为：

Y_F＝[y₀,...,y_M/2-1]

步骤2所述后参考窗谱图为：

Y_B＝[y_M/2,...,y_M-1]

其中，M表示短时傅里叶变换的时间维数，s_m,m∈[0,M-1]表示第m列STFT数据，y_m,m∈[0,M-1]表示第m列STFT谱图；

所述的s_m、y_m均为N×1维向量，其中，N为短时傅里叶变换的频率维数；

所述的S_F、S_B、Y_F、Y_B均为维矩阵；

作为优选，步骤3所述滑窗谱图为：

Y_k＝[y_kΔM,...,y_kΔM+M/2-1],k∈[1,K]

其中，M表示短时傅里叶变换的时间维数，表示相邻窗之间的时间间隔，其中/>表示向下取整操作；

y_m,m∈[0,M-1]表示第m列STFT谱图，是一个N×1维向量，N为短时傅里叶变换的频率维数，Y_k表示第k个滑窗谱图，是一个维矩阵；

作为优选，步骤4所述计算相减后前参考窗谱图与原前参考窗谱图对应元素的比值为：

步骤4所述计算相减后后参考窗谱图与原后参考窗谱图对应元素的比值为：

其中，M表示短时傅里叶变换的时间维数，N为短时傅里叶变换的频率维数，K为滑窗数目；Y_F表示前参考窗谱图，Y_B表示后参考窗谱图，Y_k表示第k个滑窗谱图；Y_F[n,m]表示前参考窗谱图在第n行第m列上的元素，Y_k[n,m]表示第k个滑窗谱图中在第n行第m列上的元素，Y_B[n,m]表示后参考窗谱图在第n行第m列上的元素；P_F,k[n,m]表示相减后前参考窗谱图与原前参考窗谱图在第n行第m列上的元素的比值；P_B,k[n,m]表示相减后后参考窗谱图与原后参考窗谱图在第n行第m列上的元素的比值；

所述的Y_F、Y_B、Y_k、P_F,k、P_B,k均为维矩阵；

步骤4所述结合设置的门限确定前参考窗中微动杂波位置为：

步骤4所述结合设置的门限确定后参考窗中微动杂波位置为：

其中，T表示设置的门限值。I_F,k表示第k个前参考窗谱图微动杂波位置标志矩阵，I_B,k表示第k个后参考窗谱图微动杂波位置标志矩阵，I_F,k[n,m]表示第k个前参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素，I_B,k[n,m]第k个后参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素。I_F表示最终前参考窗谱图微动杂波位置标志矩阵，I_B表示最终后参考窗谱图微动杂波位置标志矩阵。I_F[n,m]表示最终前参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素，I_B[n,m]表示最终后参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素。其中I_F[n,m]＝1表示前参考窗谱图在第n行第m列上存在微动杂波，I_B[n,m]＝1表示后参考窗谱图在第n行第m列上存在微动杂波；

作为优选，步骤5所述去除微动杂波后前参考窗数据为：

步骤5所述去除微动杂波后后参考窗数据为：

其中，S_F表示前参考窗数据，S_B表示后参考窗数据；S_F[n,m]表示前参考窗数据在第n行第m列上的元素，S_B[n,m]表示后参考窗数据在第n行第m列上的元素。I_F表示前参考窗谱图微动杂波位置标志矩阵，I_B表示后参考窗谱图微动杂波位置标志矩阵。I_F[n,m]表示前参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素，I_B[n,m]表示后参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素。I_F[n,m]＝1表示前参考窗谱图在第n行第m列上存在微动杂波，I_B[n,m]＝1表示后参考窗谱图在第n行第m列上存在微动杂波。表示前参考窗数据中的微动杂波分量，/>表示后参考窗数据中的微动杂波分量。/>表示前参考窗数据中的微动杂波分量在第n行第m列上的元素，/>表示后参考窗数据中的微动杂波分量在第n行第m列上的元素。所以去除微动杂波后的前参考窗数据为去除微动杂波后的后参考窗数据为/>

作为优选，步骤6所述的去除微动杂波后的前参考窗数据沿时间轴进行叠加得到的前参考窗频谱为

步骤6所述的去除微动杂波后的后参考窗数据沿时间轴进行叠加得到的后参考窗频谱为

其中，N为STFT谱图的频率维数，M为STFT谱图的时间维数。表示微动杂波去除后的前参考窗数据，/>表示微动杂波去除后的后参考窗数据。/>表示微动杂波去除后的前参考窗数据中第n行第m列上的元素，/>表示微动杂波去除后的后参考窗数据中第n行第m列上的元素。X_F表示微动杂波去除后的前参考窗频谱，X_B表示微动杂波去除后的后参考窗频谱。X_F[n]表示微动杂波去除后的前参考窗频谱中的第n个元素，X_B[n]表示微动杂波去除后的后参考窗频谱中的第n个元素。C为常数。

步骤6所述去除微动杂波后的前参考窗时域数据为x_F＝ifft(X_F)。

步骤6所述去除微动杂波后的后参考窗时域数据为x_B＝ifft(X_B)。

其中，ifft(·)表示逆傅里叶变换，x_F表示去除微动杂波后的前参考窗时域数据，x_B表示去除微动杂波后的后参考窗时域数据。x_F和x_B均为N×1维向量；

步骤6所述的去除微动杂波后的时域数据为：

其中，x_F[n]表示去除微动杂波后的前参考窗时域数据中的第n个元素，x_B[n]表示去除微动杂波后的后参考窗时域数据中的第n个元素。x^(trgt)表示去除微动杂波后的完整时域数据，x^(trgt)[n]表示去除微动杂波后的完整时域数据中的第n个元素。

与现有技术相比，本发明操作简便，计算量小，在保留目标回波的同时能够有效地去除微动杂波，极大地提高了雷达的目标探测性能及在复杂环境中的适应能力。

附图说明

图1：本发明所提方法的具体流程。

图2：实施例1中回波信号的频谱。

图3：实施例1中回波信号STFT谱图。

图4：实施例1中去除微动杂波后的STFT谱图。

图5：实施例1中去除微动杂波前后的信号频谱。

图6：实施例1中去除微动杂波前后的时域信号。

图7：实施例2中回波信号的频谱。

图8：实施例2中回波信号STFT谱图。。

图9：实施例2中去除微动杂波后的STFT谱图。

图10：实施例2中去除微动杂波前后的信号频谱。

图11：实施例2中去除微动杂波前后的时域信号。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明做进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限制本发明。

下面结合附图1至11介绍本发明的具体实施方式，具体如下：

实施例1

本实施例为仿真实施例，仿真回波信号按下式进行生成

其中P和J均设为3，表示仿真回波中包含3个目标散射点回波和3个微动散射点回波。s_p和α_j分别表示第p个目标散射点和第j个微动散射点的回波幅度，最终分别以信噪比和杂噪比的形式体现出来。3个目标回波的信噪比分别设置为[-5,20,5]dB，3个微动杂波的杂噪比分别设置为[30,25,20]dB。f_d,p表示第p,0≤p≤2个目标回波的多普勒频率，设置为[800,-300,300]Hz。f_d,j表示第j,0≤j≤2个微动杂波的平动多普勒频率，设置为[-140,-100,100]Hz。l_j表示第j,0≤j≤2个微动散射点的转动半径，设置为[13.61,27.21,5.88]m。。w_j表示第j,0≤j≤2个微动散射点的转动速度，设置为[4.25π,1.40π,7.40π]rad。θ_j表示第j,0≤j≤2个微动散射点的初相，设置为[0.7π,0,1.3π]rad。λ为信号波长，设置为0.46m。仿真信号的频谱如图2所示。

利用本发明对上述仿真信号进行处理，如图1所示，一种基于STFT谱图滑窗相消的微动杂波去除方法，具体包含以下步骤：

步骤1：将原始时域采样数据进行短时傅里叶变换，得到STFT数据及STFT谱图；

步骤1所述将原始时域采样数据进行短时傅里叶变换为：

其中，N＝4000表示短时傅里叶变换的频率维数，M＝124表示短时傅里叶变换的时间维数，S表示STFT数据，是一个4000×124维的矩阵，S[m,k]表示第m行第k列STFT数据，x表示时域采样数据序列，所述时域采样数据矩阵的长度为N＝4000，x[n]表示第n个时域采样点，w表示窗函数序列，所述窗函数的有效长度为N_w＝64，并通过补零将其拓展成长度为N＝4000的向量，它满足为常数的条件，w[n-mR]表示窗函数中的第n-mR个元素，R＝32表示短时傅里叶变换相邻窗函数之间的步进长度；

将S表示为S＝[s₀,...,s_M-1]，s_m表示第m列SFFT数据；

通过Y[m,k]＝|S[m,k]|²,m∈[0,M-1],k∈[0,N-1]计算STFT谱图；

其中，Y表示STFT谱图，Y[m,k]表示STFT谱图在第m行第k列上的元素；

将Y表示为Y＝[y₀,...,y_M-1]，y_m表示第m列STFT谱图；

步骤2：将STFT数据以时间维中心为对称轴等分为前参考窗数据、后参考窗数据，将STFT谱图以时间维中心为对称轴等分为前参考窗谱图、后参考窗谱图；

步骤2所述前参考窗数据为：

S_F＝[s₀,...,s_M/2-1]

步骤2所述后参考窗数据为：

S_B＝[s_M/2,...,s_M-1]

步骤2所述前参考窗谱图为：

Y_F＝[y₀,...,y_M/2-1]

步骤2所述后参考窗谱图为：

Y_B＝[y_M/2,...,y_M-1]

其中，M＝124表示短时傅里叶变换的时间维数，s_m,m∈[0,M-1]表示第m列STFT数据，y_m,m∈[0,M-1]表示第m列STFT谱图；

所述的s_m、y_m均为N×1维向量，其中，N＝4000为短时傅里叶变换的频率维数；

所述的S_F、S_B、Y_F、Y_B均为4000×62维矩阵；

步骤3：设置滑窗数目K＝6，从STFT谱图中截取得到K个滑窗谱图；

步骤3所述滑窗谱图为：

Y_k＝[y_kΔM,...,y_kΔM+M/2-1],k∈[1,K]

其中，M＝124表示短时傅里叶变换的时间维数，表示相邻窗之间的时间间隔，其中/>表示向下取整操作；

y_m,m∈[0,M-1]表示第m列STFT谱图，是一个N×1维向量，N＝4000为短时傅里叶变换的频率维数，Y_k表示第k个滑窗谱图，是一个4000×62维矩阵；

步骤4：将前参考窗谱图与各滑窗谱图进行对应元素相减得到相减后前参考窗谱图，计算相减后前参考窗谱图与原前参考窗谱图对应元素的比值；将后参考窗谱图与各滑窗谱图进行对应元素相减得到相减后后参考窗谱图，计算相减后后参考窗谱图与原后参考窗谱图对应元素的比值；结合设置的门限确定前参考窗谱图中微动杂波位置；结合设置的门限确定后参考窗谱图中微动杂波位置；

步骤4所述计算相减后前参考窗谱图与原前参考窗谱图对应元素的比值为：

步骤4所述计算相减后后参考窗谱图与原后参考窗谱图对应元素的比值为：

其中，M＝124表示短时傅里叶变换的时间维数，N＝4000为短时傅里叶变换的频率维数，K＝6为滑窗数目；Y_F表示前参考窗谱图，Y_B表示后参考窗谱图，Y_k表示第k个滑窗谱图；Y_F[n,m]表示前参考窗谱图在第n行第m列上的元素，Y_k[n,m]表示第k个滑窗谱图中在第n行第m列上的元素，Y_B[n,m]表示后参考窗谱图在第n行第m列上的元素；P_F,k[n,m]表示相减后前参考窗谱图与原前参考窗谱图在第n行第m列上的元素的比值；P_B,k[n,m]表示相减后后参考窗谱图与原后参考窗谱图在第n行第m列上的元素的比值；

所述的Y_F、Y_B、Y_k、P_F,k、P_B,k均为4000×62维矩阵；

步骤4所述结合设置的门限确定前参考窗中微动杂波位置为：

步骤4所述结合设置的门限确定后参考窗中微动杂波位置为：

其中，T＝0.65表示设置的门限值。I_F,k表示第k个前参考窗谱图微动杂波位置标志矩阵，I_B,k表示第k个后参考窗谱图微动杂波位置标志矩阵，I_F,k[n,m]表示第k个前参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素，I_B,k[n,m]第k个后参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素。I_F表示最终前参考窗谱图微动杂波位置标志矩阵，I_B表示最终后参考窗谱图微动杂波位置标志矩阵。I_F[n,m]表示最终前参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素，I_B[n,m]表示最终后参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素。其中I_F[n,m]＝1表示前参考窗谱图在第n行第m列上存在微动杂波，I_B[n,m]＝1表示后参考窗谱图在第n行第m列上存在微动杂波；

步骤5：将微动杂波所在位置上的前参考窗数据置零，得到去除微动杂波后的前参考窗数据；将微动杂波所在位置上的后参考窗数据置零，得到去除微动杂波后的后参考窗数据；

步骤5所述去除微动杂波后前参考窗数据为：

步骤5所述去除微动杂波后后参考窗数据为：

其中，S_F表示前参考窗数据，S_B表示后参考窗数据；S_F[n,m]表示前参考窗数据在第n行第m列上的元素，S_B[n,m]表示后参考窗数据在第n行第m列上的元素。I_F表示前参考窗谱图微动杂波位置标志矩阵，I_B表示后参考窗谱图微动杂波位置标志矩阵。I_F[n,m]表示前参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素，I_B[n,m]表示后参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素。I_F[n,m]＝1表示前参考窗谱图在第n行第m列上存在微动杂波，I_B[n,m]＝1表示后参考窗谱图在第n行第m列上存在微动杂波。表示前参考窗数据中的微动杂波分量，/>表示后参考窗数据中的微动杂波分量。/>表示前参考窗数据中的微动杂波分量在第n行第m列上的元素，/>表示后参考窗数据中的微动杂波分量在第n行第m列上的元素。所以去除微动杂波后的前参考窗数据为去除微动杂波后的后参考窗数据为/>

步骤6：将去除微动杂波后的前参考窗数据沿时间轴进行叠加得到前参考窗频谱，对前参考窗频谱进行逆傅里叶变换得到去除微动杂波后的前参考窗时域数据；将去除微动杂波后的后参考窗数据沿时间轴进行叠加得到后参考窗频谱，对后参考窗频谱进行逆傅里叶变换得到去除微动杂波后的后参考窗时域数据；将去除微动杂波后的前参考窗时域数据和去除微动杂波后的后参考窗时域数据进行组合得到去除微动杂波后的完整时域数据。

步骤6所述的去除微动杂波后的前参考窗数据沿时间轴进行叠加得到的前参考窗频谱为

步骤6所述的去除微动杂波后的后参考窗数据沿时间轴进行叠加得到的后参考窗频谱为

其中，N＝4000为STFT谱图的频率维数，M＝124为STFT谱图的时间维数。表示微动杂波去除后的前参考窗数据，/>表示微动杂波去除后的后参考窗数据。/>表示微动杂波去除后的前参考窗数据中第n行第m列上的元素，/>表示微动杂波去除后的后参考窗数据中第n行第m列上的元素。X_F表示微动杂波去除后的前参考窗频谱，X_B表示微动杂波去除后的后参考窗频谱。X_F[n]表示微动杂波去除后的前参考窗频谱中的第n个元素，X_B[n]表示微动杂波去除后的后参考窗频谱中的第n个元素。C＝1.08为常数。

步骤6所述去除微动杂波后的前参考窗时域数据为x_F＝ifft(X_F)。

步骤6所述去除微动杂波后的后参考窗时域数据为x_B＝ifft(X_B)。

其中，ifft(·)表示逆傅里叶变换，x_F表示去除微动杂波后的前参考窗时域数据，x_B表示去除微动杂波后的后参考窗时域数据。x_F和x_B均为4000×1维向量；

步骤6所述的去除微动杂波后的时域数据为：

其中，x_F[n]表示去除微动杂波后的前参考窗时域数据中的第n个元素，x_B[n]表示去除微动杂波后的后参考窗时域数据中的第n个元素。x^(trgt)表示去除微动杂波后的完整时域数据，x^(trgt)[n]表示去除微动杂波后的完整时域数据中的第n个元素。

实施例2

本实施例通过外场实验来验证本发明的有效性。本实施例利用外辐射源雷达系统在风电场区内对大疆无人机进行监测，因此雷达回波中将存在由风电机组叶片引入的微动杂波。所用外辐射源雷达系统是以地面数字多媒体广播信号(DTMB)作为机会照射源，信号中心频率为722MHz，带宽7.56MHz。慢时间维采样率为1800Hz，风电机组叶片长度约为56.8m，转速约为13.32rpm。处理所用积累时间为4.44s。如图7所示为原始回波信号，可以看出无人机多普勒频率约为-60Hz，落在微动杂波区内。利用所提算法对该回波数据进行处理，由于本实施例和实施例1为同一发明技术构思，在本实施例中未详尽描述的内容，请参见实施例1。本实施例具体包括以下步骤：

步骤1：将原始时域采样数据进行短时傅里叶变换，得到STFT数据及STFT谱图。在本实施例中，所用窗函数是有效长度N_w＝256的汉明窗。相邻窗之间的步进值R＝128。STFT数据和STFT谱图的频率维数均为N＝8000，时间维数均为M＝60。如图8所示为步骤1得到的回波信号STFT谱图。

步骤2：以时间维中心为对称轴将STFT数据和STFT谱图分别等分为前、后两段，称为前、后参考窗STFT数据及前、后参考窗STFT谱图。在本实施例中前、后参考窗STFT数据及前、后参考窗STFT谱图均为8000×30的矩阵。

步骤3：设置滑窗数目，并从STFT谱图中截取得到对应数目的滑窗STFT谱图。在本实施例中设置滑窗数目K＝6，则相邻窗函数之间的时间间隔ΔM＝4，截取的滑窗STFT谱图为均为8000×30的矩阵。

步骤4：将前、后参考窗谱图分别与各滑窗STFT谱图进行相减，并根据相减前后谱图强度的变化情况确定前、后参考窗中微动杂波位置。在本实施例中设置相减前后强度变化的门限T＝0.7。

步骤5：将微动杂波所在位置上的前、后参考窗STFT数据置零，得到去除微动杂波后的前、后参考窗STFT数据。如图9所示为微动去除后的STFT谱图。

步骤6：将去除微动杂波后的前、后参考窗STFT数据分别沿时间轴进行叠加，然后对叠加后的结果分别进行逆傅里叶变换，并将两个逆傅里叶变换的结果进行组合得到去除微动杂波后的时域数据。图10、图11分别为微动杂波去除前后的频谱和时域回波图。从中可以看出风车引起的微动杂波被有效地抑制，同时无人机回波得到明显凸显，证明了所提算法的有效性。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (7)

1.一种基于STFT谱图滑窗相消的微动杂波去除方法，其特征在于，包含以下步骤：

步骤1：将原始时域采样数据进行短时傅里叶变换，得到STFT数据及STFT谱图；

步骤2：将STFT数据以时间维中心为对称轴等分为前参考窗数据、后参考窗数据，将STFT谱图以时间维中心为对称轴等分为前参考窗谱图、后参考窗谱图；

步骤3：设置滑窗数目K，从STFT谱图中截取得到K个滑窗谱图；

步骤4：将前参考窗谱图与各滑窗谱图进行对应元素相减得到相减后前参考窗谱图，计算相减后前参考窗谱图与原前参考窗谱图对应元素的比值；将后参考窗谱图与各滑窗谱图进行对应元素相减得到相减后后参考窗谱图，计算相减后后参考窗谱图与原后参考窗谱图对应元素的比值；结合设置的门限确定前参考窗谱图中微动杂波位置；结合设置的门限确定后参考窗谱图中微动杂波位置；

步骤5：将微动杂波所在位置上的前参考窗数据置零，得到去除微动杂波后的前参考窗数据；将微动杂波所在位置上的后参考窗数据置零，得到去除微动杂波后的后参考窗数据；

步骤6：将去除微动杂波后的前参考窗数据沿时间轴进行叠加得到前参考窗频谱，对前参考窗频谱进行逆傅里叶变换得到去除微动杂波后的前参考窗时域数据；将去除微动杂波后的后参考窗数据沿时间轴进行叠加得到后参考窗频谱，对后参考窗频谱进行逆傅里叶变换得到去除微动杂波后的后参考窗时域数据；将去除微动杂波后的前参考窗时域数据和去除微动杂波后的后参考窗时域数据进行组合得到去除微动杂波后的完整时域数据。

2.根据权利要求1所述的基于STFT谱图滑窗相消的微动杂波去除方法，其特征在于，步骤1所述将原始时域采样数据进行短时傅里叶变换为：

其中，N表示短时傅里叶变换的频率维数，M表示短时傅里叶变换的时间维数，S表示STFT数据，是一个N×M维的矩阵，S[m,k]表示第m行第k列STFT数据，x表示时域采样数据序列，所述时域采样数据矩阵的长度为N，x[n]表示第n个时域采样点，w表示窗函数序列，所述窗函数的有效长度为N_w，并通过补零将其拓展成长度为N的向量，它满足为常数的条件，w[n-mR]表示窗函数中的第n-mR个元素，R表示短时傅里叶变换相邻窗函数之间的步进长度；

将S表示为S＝[s₀,...,s_M-1]，s_m表示第m列SFFT数据；

通过Y[m,k]＝|S[m,k]|²,m∈[0,M-1],k∈[0,N-1]计算STFT谱图；

其中，Y表示STFT谱图，Y[m,k]表示STFT谱图在第m行第k列上的元素；

将Y表示为Y＝[y₀,...,y_M-1]，y_m表示第m列STFT谱图。

3.根据权利要求1所述的基于STFT谱图滑窗相消的微动杂波去除方法，其特征在于，步骤2所述前参考窗数据为：

S_F＝[s₀,...,s_M/2-1]

步骤2所述后参考窗数据为：

S_B＝[s_M/2,...,s_M-1]

步骤2所述前参考窗谱图为：

Y_F＝[y₀,...,y_M/2-1]

步骤2所述后参考窗谱图为：

Y_B＝[y_M/2,...,y_M-1]

其中，M表示短时傅里叶变换的时间维数，s_m,m∈[0,M-1]表示第m列STFT数据，y_m,m∈[0,M-1]表示第m列STFT谱图；

所述的s_m、y_m均为N×1维向量，其中，N为短时傅里叶变换的频率维数；

所述的S_F、S_B、Y_F、Y_B均为维矩阵。

4.根据权利要求1所述的基于STFT谱图滑窗相消的微动杂波去除方法，其特征在于，步骤3所述滑窗谱图为：

Y_k＝[y_kΔM,...,y_kΔM+M/2-1],k∈[1,K]

其中，M表示短时傅里叶变换的时间维数，ΔM＝M/[2(K+1)]表示相邻窗之间的时间间隔，其中·表示向下取整操作；

y_m,m∈[0,M-1]表示第m列STFT谱图，是一个N×1维向量，N为短时傅里叶变换的频率维数，Y_k表示第k个滑窗谱图，是一个维矩阵。

5.根据权利要求1所述的基于STFT谱图滑窗相消的微动杂波去除方法，其特征在于，步骤4所述计算相减后前参考窗谱图与原前参考窗谱图对应元素的比值为：

步骤4所述计算相减后后参考窗谱图与原后参考窗谱图对应元素的比值为：

其中，M表示短时傅里叶变换的时间维数，N为短时傅里叶变换的频率维数，K为滑窗数目；Y_F表示前参考窗谱图，Y_B表示后参考窗谱图，Y_k表示第k个滑窗谱图；Y_F[n,m]表示前参考窗谱图在第n行第m列上的元素，Y_k[n,m]表示第k个滑窗谱图中在第n行第m列上的元素，Y_B[n,m]表示后参考窗谱图在第n行第m列上的元素；P_F,k[n,m]表示相减后前参考窗谱图与原前参考窗谱图在第n行第m列上的元素的比值；P_B,k[n,m]表示相减后后参考窗谱图与原后参考窗谱图在第n行第m列上的元素的比值；

所述的Y_F、Y_B、Y_k、P_F,k、P_B,k均为维矩阵；

步骤4所述结合设置的门限确定前参考窗中微动杂波位置为：

步骤4所述结合设置的门限确定后参考窗中微动杂波位置为：

其中，I_F,k表示第k个前参考窗谱图微动杂波位置标志矩阵，I_B,k表示第k个后参考窗谱图微动杂波位置标志矩阵，I_F,k[n,m]表示第k个前参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素，I_B,k[n,m]第k个后参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素；I_F表示最终前参考窗谱图微动杂波位置标志矩阵，I_B表示最终后参考窗谱图微动杂波位置标志矩阵；I_F[n,m]表示最终前参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素，I_B[n,m]表示最终后参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素；其中I_F[n,m]＝1表示前参考窗谱图在第n行第m列上存在微动杂波，I_B[n,m]＝1表示后参考窗谱图在第n行第m列上存在微动杂波。

6.根据权利要求1所述的基于STFT谱图滑窗相消的微动杂波去除方法，其特征在于，步骤5所述去除微动杂波后前参考窗数据为：

步骤5所述去除微动杂波后后参考窗数据为：

其中，S_F表示前参考窗数据，S_B表示后参考窗数据；S_F[n,m]表示前参考窗数据在第n行第m列上的元素，S_B[n,m]表示后参考窗数据在第n行第m列上的元素；I_F表示前参考窗谱图微动杂波位置标志矩阵，I_B表示后参考窗谱图微动杂波位置标志矩阵；I_F[n,m]表示前参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素，I_B[n,m]表示后参考窗谱图微动杂波位置标志矩阵在第n行第m列上的元素；I_F[n,m]＝1表示前参考窗谱图在第n行第m列上存在微动杂波，I_B[n,m]＝1表示后参考窗谱图在第n行第m列上存在微动杂波；表示前参考窗数据中的微动杂波分量，/>表示后参考窗数据中的微动杂波分量；/>表示前参考窗数据中的微动杂波分量在第n行第m列上的元素，/>表示后参考窗数据中的微动杂波分量在第n行第m列上的元素；所以去除微动杂波后的前参考窗数据为去除微动杂波后的后参考窗数据为/>

7.根据权利要求1所述的基于STFT谱图滑窗相消的微动杂波去除方法，其特征在于，步骤6所述的去除微动杂波后的前参考窗数据沿时间轴进行叠加得到的前参考窗频谱为：

步骤6所述的去除微动杂波后的后参考窗数据沿时间轴进行叠加得到的后参考窗频谱为

其中，N为STFT谱图的频率维数，M为STFT谱图的时间维数；表示微动杂波去除后的前参考窗数据，/>表示微动杂波去除后的后参考窗数据；/>表示微动杂波去除后的前参考窗数据中第n行第m列上的元素，/>表示微动杂波去除后的后参考窗数据中第n行第m列上的元素；X_F表示微动杂波去除后的前参考窗频谱，X_B表示微动杂波去除后的后参考窗频谱；X_F[n]表示微动杂波去除后的前参考窗频谱中的第n个元素，X_B[n]表示微动杂波去除后的后参考窗频谱中的第n个元素；C为常数；

步骤6所述去除微动杂波后的前参考窗时域数据为x_F＝ifft(X_F)；

步骤6所述去除微动杂波后的后参考窗时域数据为x_B＝ifft(X_B)；

其中，ifft(·)表示逆傅里叶变换，x_F表示去除微动杂波后的前参考窗时域数据，x_B表示去除微动杂波后的后参考窗时域数据；x_F和x_B均为N×1维向量；

步骤6所述的去除微动杂波后的时域数据为：

其中，x_F[n]表示去除微动杂波后的前参考窗时域数据中的第n个元素，x_B[n]表示去除微动杂波后的后参考窗时域数据中的第n个元素；x^(trgt)表示去除微动杂波后的完整时域数据，x^(trgt)[n]表示去除微动杂波后的完整时域数据中的第n个元素。