CN113238210B - 基于微多普勒效应的悬停无人机特征检测方法 - Google Patents

基于微多普勒效应的悬停无人机特征检测方法 Download PDF

Info

Publication number
CN113238210B
CN113238210B CN202110481764.7A CN202110481764A CN113238210B CN 113238210 B CN113238210 B CN 113238210B CN 202110481764 A CN202110481764 A CN 202110481764A CN 113238210 B CN113238210 B CN 113238210B
Authority
CN
China
Prior art keywords
test data
echo
time
frequency
training data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110481764.7A
Other languages
English (en)
Other versions
CN113238210A (zh
Inventor
李明
王皓
左磊
李祥朋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN202110481764.7A priority Critical patent/CN113238210B/zh
Publication of CN113238210A publication Critical patent/CN113238210A/zh
Application granted granted Critical
Publication of CN113238210B publication Critical patent/CN113238210B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/02Systems using reflection of radio waves, e.g. primary radar systems; Analogous systems
    • G01S13/04Systems determining presence of a target
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种基于微多普勒效应的悬停无人机特征检测方法,主要解决现有技术在低重频、低信杂比下检测效果差的问题。其方案是:接收悬停无人机的雷达回波,并对其进行下变频和脉冲压缩处理;将脉冲压缩信号每一距离单元的慢时间维数据分为训练集和测试集;分别计算训练数据和测试数据的时域Hurst指数、频谱熵值、时频脊的均值、时频脊的方差和倒谱熵值;用训练数据的特征值组成训练特征向量,训练支撑向量机,生成支撑向量机检测器;用测试数据的特征值组成测试特征向量,输入到支撑向量机检测器,对悬停无人机进行特征检测。本发明利用金属转子微多普勒特征提高了对悬停无人机的特征检测效果,可用于低重频、强地杂波环境下的目标检测。

Description

基于微多普勒效应的悬停无人机特征检测方法
技术领域
本发明属于雷达技术领域,更进一步涉及一种悬停无人机特征检测方法。本发明可用于低重频、强地杂波环境下的目标检测。
背景技术
近年来,随着无人机技术的快速发展,无人机的管控面临着巨大的困难。无人机这种低慢小目标由于其飞行高度“低”、飞行速度“慢”、雷达反射截面积“小”,很容易被淹没在强烈的地杂波中,传统的频域滤波检测方法检测性能严重下降,难以完成对无人机等低慢小目标的检测。
北京航空航天大学在其申请的专利文献“一种基于双发双收相干雷达的无人机旋翼长度和转速的估计方法”(申请号201810779125.7,公开号CN 108957443 A)中公开了一种对无人机的旋翼长度和转速进行精确估计的方法。该方法首先通过对发射信号进行相位补偿,形成波束对准目标;然后,利用两路接收天线收集旋翼回波信号,并将两路回波信号输入相干器,得到相干信号,通过相干信号的自相关函数提取旋翼的旋转半周期;最后,利用时频谱特定时刻具有幅度极大值的频率的差值,估计无人机旋翼的长度和转速,实现无人机的特征检测和识别。但是,由于无人机旋翼回波的主瓣极窄,当雷达的脉冲重复频率较小时,使用该方法会产生雷达采不到无人机旋翼回波主瓣的问题,导致其无法适用于低重频情况下的无人机参数估计和目标检测。
方鑫在其发表的论文“小型无人机目标雷达探测关键技术研究”(电子科技大学2019博士论文)中提出了一种基于旋翼微多普勒特征的无人机目标检测方法。该方法根据不同叶片微多普勒频率差值与雷达载波频率、叶片长度、叶片旋转速度以及叶片旋转初始角之间的关系,给出了根据微多普勒特征区分不同叶片回波信号的条件,基于时频分析、HRT和二维CFAR实现了旋翼叶片的微多普勒特征提取和无人机目标的检测。该方法由于没有考虑杂波的影响,且使用旋翼微多普勒在时频域的正弦曲线对无人机进行检测,而实际数据中无人机旋翼在时频域的曲线仅在高信噪比、高脉冲采样频率的情况下才能被观察到的,因而导致其在低重频、低信杂比情况下的检测效果较差。
发明内容
本发明的目的在于针对上述已有技术的不足,提出一种基于微多普勒效应的悬停无人机特征检测方法,以提取无人机金属转子和旋翼的微多普勒特征,在低重频、低信杂比情况下实现对无人机的有效检测。
为实现上述目的,本发明的技术方案包括如下:
(1)对悬停无人机的雷达回波信号在快时间域进行下变频和脉冲压缩处理,并提取每一个距离单元的慢时间域信号;
(2)将慢时间域数据分为目标单元回波和非目标单元回波两部分,取一半目标单元回波和一半非目标单元回波组成训练数据集,将剩余的回波数据作为测试数据集;
(3)根据训练数据集和测试数据集的慢时间域数据,在不同变换域中计算表征无人机金属转子和旋翼微多普勒特性的特征值:
(3a)在时域计算训练数据的时域赫斯特指数H0和测试数据的时域赫斯特指数H1
(3b)在频域计算训练数据的频谱熵值Ex,0和测试数据的频谱熵值Ex,1
(3c)在时频域计算训练数据时频脊的均值mtfr,0和方差VARtfr,0,同时计算测试数据时频脊的均值mtfr,1和方差VARtfr,1
(3d)在倒频域计算训练数据的倒谱熵值Eq,0和测试数据的倒谱熵值Eq,1
(4)用上述计算出训练数据集中的所有特征值组成训练特征向量,对支撑向量机进行训练,生成支撑向量机检测器;
(5)使用上述计算出测试数据集的所有特征值组成特征向量,输入到支撑向量机检测器,判断是否有目标存在,实现对悬停无人机的特征检测。
本发明利用时域赫斯特指数、频谱熵值、时频脊的均值和方差、倒谱熵值这5个微多普勒特征对悬停无人机进行特征检测,相比现有技术具有如下优点:
第一,由于这些特征既包含旋翼的微多普勒效应也包含金属转子的微多普勒效应,因此,目标单元回波的特征和非目标单元回波的特征在低重频条件下依旧有明显的区别,可在低重频条件下对悬停无人机进行有效的检测,克服了现有基于旋翼微多普勒的无人机目标检测技术在低重频条件下检测效果差的问题。
第二,由于微多普勒特征受杂波的影响较小,因此,微多普勒特征能够有效区分纯杂波信号和淹没在杂波中的无人机信号,在信杂比较低时也能对悬停无人机进行有效的检测,克服了现有技术难以对淹没在杂波中的无人机进行有效检测的问题,提高了在强杂波环境下的检测效果。
附图说明
图1为本发明的实现流程图;
图2为仿真的目标单元回波;
图3为用本发明对图2进行特征检测的仿真实验结果图。
具体实施方式
下面结合附图,对本发明的实施例和效果作进一步详细描述。
参照图1,对实施的实现步骤如下:
步骤1,对悬停无人机的雷达回波信号在快时间域进行下变频和脉冲压缩处理。
(1.1)获取悬停无人机的雷达回波。
无人机为一悬停的四旋翼无人机,雷达位于一处高楼的楼顶,雷达视线方向正对着无人机进行观测,从雷达的接收机获取悬停无人机的回波信号sr
sr=s0+sm+sR+sc+sn
其中,s0为无人机机身的回波,sm为无人机金属转子的回波,sR为无人机旋翼的回波,sc为杂波,sn为噪声;
(1.2)获取悬停无人机的下变频信号:
在快时间域内对悬停无人机的雷达回波信号进行下变频,得到悬停无人机的下变频信号h:
h=sr exp(-j2πztr),
其中,exp表示以自然常数e为底的指数操作,j表示虚数单位符号,z表示根据工程需求设定的雷达机动目标回波信号的载频,tr表示雷达接收回波信号快时间域的采样时间点,tr的取值范围为
Figure BDA0003049521970000031
fs表示根据工程需求设定的雷达接收机动回波信号的采样频率;
(1.3)获取悬停无人机的脉压信号:
在快时间域内对下变频信号在距离维进行脉冲压缩,得到悬停无人机的脉冲压缩信号g:
Figure BDA0003049521970000032
其中,
Figure BDA0003049521970000041
表示卷积操作,rect(·)表示矩形窗因子,当
Figure BDA0003049521970000042
时矩形窗因子取值为0,当
Figure BDA0003049521970000043
时矩形窗因子取值为1,Tp表示雷达接收的回波信号的脉冲宽度,|·|表示取绝对值操作,k表示雷达接收的回波信号的调频斜率。
步骤2,提取每一个距离单元的慢时间域信号。
在悬停无人机的脉压信号中,提取每一个距离单元的慢时间域信号,包括非目标单元和目标单元的慢时间域回波信号,分别表示如下:
非目标单元的慢时间域回波信号x(n)为:
x(n)=sc+sn
其中,n∈[1,N],N为单帧数据的脉冲个数,sc为杂波,sn为噪声;
目标单元的慢时间域回波信号x(n)为:
Figure BDA0003049521970000044
其中,a和b为数值不同的两个常数,t=nTr为慢时间域采样点,Tr为脉冲重复时间,s0为机身回波,sR,k(t)为第k个旋翼的回波,sm,k(t)为第k个金属转子的回波。
所述的无人机旋翼回波sR,k(t)和金属转子的回波的公式分别表示如下:
Figure BDA0003049521970000045
Figure BDA0003049521970000046
其中,L为旋翼的长度,fd为无人机径向运动的多普勒频率,λ为雷达信号的波长,R0为悬停无人机到雷达的初始距离,nb为单个旋翼的叶片数量,β为无人机相对雷达视线方向的俯仰角,ψt,k=2πfΩk为旋翼的旋转角,即雷达与旋翼所在的直线和旋翼叶片所在直线之间的夹角,fΩ为旋转速度,θk为第k个旋翼的初始旋转角,d为常数,J0(·)为第一类0阶贝塞尔函数,r为金属转子的半径,γt,k为金属转子表面的某散射点P所在切线和雷达与旋翼所在直线的夹角,ψt,k为金属转子表面的散射点P的旋转角,U1(t)和U2(t)为周期矩形脉冲信号:
Figure BDA0003049521970000051
U2(t)=1-U1(t)。
步骤3,对慢时间域数据进行分集。
将慢时间域数据分为目标单元回波和非目标单元回波两部分,取一半目标单元回波和一半非目标单元回波组成训练数据集,将剩余的回波数据作为测试数据集。
步骤4,根据训练数据集和测试数据集的慢时间域数据,在不同变换域中计算表征无人机金属转子和旋翼微多普勒特性的特征值。
表征无人机金属转子和旋翼微多普勒特性的特征值,包括时域赫斯特指数H、频谱熵值Ef、时频脊的均值mtfr、时频脊的方差VARtfr和倒谱熵值Eq。这些特征值可以表征悬停无人机在不同变换域的微多普勒特征。
由于地面散射体都是固定不动的,而悬停无人机的旋翼和金属转子都是在快速转动的,所以地杂波的波动性较弱,而悬停无人机的回波的波动性较强。因此,目标单元回波与非目标单元回波的时域波动性有明显的差异,使用时域赫斯特指数H可以计量时域信号的波动性,作为悬停无人机特征检测的一个特征。
由于目标单元回波中有悬停无人机的微多普勒分量,所以其频谱的能量分布较为分散,频谱的混乱度较大。而非目标单元回波中没有微多普勒分量,其频谱的能量主要集中于杂波所在的零频处,混乱度较小。因此,目标单元和非目标单元的频谱混乱度有明显的差异,使用频谱熵值Ef可以计量信号频谱的混乱度,作为悬停无人机特征检测的一个特征。
由于目标单元回波中有悬停无人机的微多普勒分量,所以当短时傅里叶的窗长较长时,目标单元回波的时频谱中会出现关于零频对称的相互平行的微多普勒谱线,即金属转子的微多普勒时频脊。对于目标单元回波,除零频外,每一时刻瞬时频率的幅度最大值都在上述时频脊上。对于非目标单元回波,除零频外,每一时刻瞬时频率的幅度最大值都是随机分布的噪声分量。因此,可以提取每一时刻瞬时频率的幅度最大值所在的频率,计算其均值和方差,即时频脊的均值mtfr和时频脊的方差VARtfr,作为悬停无人机特征检测的特征。
由于目标单元回波的频谱中有周期性出现的微多普勒谱线,所以其倒谱的能量分布较为分散,倒谱的混乱度较大。而非目标单元回波中没有微多普勒分量,其倒谱的能量主要集中在倒频为零处,混乱度较小。因此,目标单元和非目标单元的倒谱混乱度有明显的差别,使用倒谱熵值Eq可以计量信号倒谱的混乱度,作为悬停无人机特征检测的一个特征。
上述特征值的计算分别如下:
(4.1)在时域计算训练数据的时域赫斯特指数H0和测试数据的时域赫斯特指数H1
(4.1.1)分别从训练数据集中取慢时间域信号x0(n),计算训练数据的累积离差Y0(n),从测试数据集中取慢时间域信号x1(n),计算测试数据的累积离差Y1(n):
Figure BDA0003049521970000061
Figure BDA0003049521970000062
其中,n∈[1,N],N为单帧数据的脉冲个数,<·>表示取均值的操作;
(4.1.2)分别从训练数据的累积离差序列Y0(n)的开头和末尾开始划分出
Figure BDA0003049521970000063
个连续且长度为m的短序列,构成2Nm个短序列,同时,对测试数据的累积离差序列Y1(n)执行同样的划分操作,也构成2Nm个短序列;
(4.1.3)使用最小二乘法分别对序列Y0(n)和序列Y1(n)划分出的短序列进行拟合,分别得到训练数据的拟合多项式y0(i)和测试数据的拟合多项式y1(i),其中,i的取值范围为[1,2Nm];
(4.1.4)分别对训练数据和测试数据进行消除趋势估计,获得训练数据的拟合方差Yυ,0(s,m)和测试数据的拟合方差Yυ,1(s,m):
Figure BDA0003049521970000071
Figure BDA0003049521970000072
其中,s为拟合多项式的序号,s∈[1,2Nm];
(4.1.5)对(4.1.4)得到的两个拟合方差分别取均值,再开方,得到训练数据的二阶波动函数F0(m)和测试数据的二阶波动函数F1(m):
Figure BDA0003049521970000073
Figure BDA0003049521970000074
(4.1.6)根据训练数据的二阶波动函数F0(m)和测试数据的二阶波动函数F1(m),使用最小二乘法,分别计算训练数据的时域赫斯特指数H0和测试数据的时域赫斯特指数H1
Figure BDA0003049521970000075
Figure BDA0003049521970000076
(4.2)在频域计算训练数据的频谱熵值Ef,0和测试数据的频谱熵值Ef,1
(4.2.1)利用傅里叶变换公式,分别计算训练数据x0(n)的频谱X0(f)和测试数据x1(n)的频谱X1(f):
Figure BDA0003049521970000077
Figure BDA0003049521970000081
其中,f为频率,f∈[-fr/2,fr/2],fr为雷达的脉冲重复频率;
(4.2.2)分别对训练数据的频谱X0(f)和测试数据的频谱X1(f)进行归一化,得到训练数据的归一化频谱X0,1(f)和测试数据的归一化频谱X1,1(f):
Figure BDA0003049521970000082
Figure BDA0003049521970000083
(4.2.3)分别计算训练数据的频谱熵值Ef,0和测试数据的频谱熵值Ef,1
Figure BDA0003049521970000084
Figure BDA0003049521970000085
(4.3)在时频域计算训练数据时频脊的均值mtfr,0和方差VARtfr,0,同时计算测试数据时频脊的均值mtfr,1和方差VARtfr,1
(4.3.1)将短时傅里叶变换的短时窗设置为汉明窗,为使金属转子的微多普勒特征易于观察和提取,设窗长应大于旋翼旋转周期的2倍,即:
Lw>2T,
其中,Lw为短时傅里叶变换的窗长,T为无人机旋翼和金属转子的旋转周期;
(4.3.2)分别对训练数据x0(n)和测试数据x1(n)作短时傅里叶变换,得到训练数据的时频谱STFTx0(n,f)和测试数据的时频谱STFTx1(n,f):
Figure BDA0003049521970000086
Figure BDA0003049521970000087
其中,n∈[1,N],f∈[-fr/2,fr/2],w(·)为短时窗函数;
(4.3.3)分别对训练数据和测试数据的时频谱进行截取,得到去除机身信号的训练数据时频谱STFTx0,0(n,f)和去除机身信号的测试数据时频谱STFTx1,0(n,f):
Figure BDA0003049521970000091
Figure BDA0003049521970000092
(4.3.4)根据去除机身信号的时频谱,分别计算训练数据中无人机金属转子的时频脊R0(n)和测试数据中无人机金属转子的时频脊R1(n):
R0(n)={|fm||STFTx0,0(n,fm)=max[STFTx0,0(n,f)]},
R1(n)={|fm||STFTx1,0(n,fm)=max[STFTx1,0(n,f)]},
(4.3.5)根据(4.3.4)的结果,分别计算训练数据时频脊的均值mtfr,0和测试数据时频脊的均值mtfr,1
Figure BDA0003049521970000093
Figure BDA0003049521970000094
(4.3.6)根据(4.3.4)和(4.3.5)的结果,分别计算训练数据时频脊的方差VARtfr,0和测试数据时频脊的方差VARtfr,1
Figure BDA0003049521970000095
Figure BDA0003049521970000096
(4.4)在倒频域计算训练数据的倒谱熵值Eq,0和测试数据的倒谱熵值Eq,1
(4.4.1)利用倒谱变换公式,分别计算训练数据x0(n)的倒谱C0(q)和测试数据x1(n)的倒谱C1(q):
C0(q)=|IFFT{log|FFT[x0(n)]|}|,
C1(q)=|IFFT{log|FFT[x1(n)]|}|,
其中,q为倒频率,FFT(·)为傅里叶变换,IFFT(·)为逆傅里叶变换;
(4.4.2)分别对训练数据和测试数据的倒谱进行截取,得到去除机身信号的训练数据倒谱C0,0(q)和去除机身信号的测试数据倒谱C1,0(q):
C0,0(q)=C0(q),q∈(T/3,NTr/2),
C1,0(q)=C1(q),q∈(T/3,NTr/2);
(4.4.3)对截取的倒谱进行归一化,分别计算训练数据的归一化倒谱C0,1(q)和测试数据的归一化倒谱C1,1(q):
Figure BDA0003049521970000101
Figure BDA0003049521970000102
(4.4.4)根据(4.4.3)的结果,分别计算训练数据的倒谱熵值Eq,0和测试数据的倒谱熵值Eq,1
Figure BDA0003049521970000103
Figure BDA0003049521970000104
步骤5,根据训练数据集中的所有特征值,对支撑向量机进行训练,生成支撑向量机检测器。
(5.1)分别使用训练数据集中每一个慢时间域信号的特征值,组成训练特征向量f0,i
f0,i=(H0,Ef,0,mtfr,0,VARtfr,0,Eq,0),
其中,i∈[1,M0],M0为训练数据集中慢时间域回波的数量;
(5.2)为训练特征向量f0,i标注类别,得到训练样本集Ω0
Ω0={(f0,i,yi)|i=1,2,...,M0},
其中,yi为慢时间域信号的类别,若慢时间域信号来自目标单元,则yi=1,若慢时间域信号来自非目标单元,则yi=-1;
(5.3)根据训练样本集构造一个最优超平面g(f)=ωTf+b作为决策曲面,使得目标单元回波和非目标单元回波样本之间的隔离边缘最大化,该最优超平面需满足如下条件:
yiTf0,i+b)≥1,i=1,2,...,M0
其中,ω为超平面的法向量,b为超平面的位移项,(·)T为转置操作;
(5.4)根据(5.3)中的限制条件,将最优超平面参数的求解问题转化为如下的二次规划问题:
Figure BDA0003049521970000111
再通过拉格朗日乘数法,将上述二次规划问题转化为对偶问题Q(α):
Figure BDA0003049521970000112
其中,
Figure BDA0003049521970000113
为拉格朗日系数,αi为拉格朗日系数α中的第i个值,αi≥0,i∈[1,M0];
(5.5)对(5.4)中的对偶问题Q(α)进行求解,得到拉格朗日系数α的解为
Figure BDA0003049521970000114
再由拉格朗日系数的解α*计算超平面法向量的解ω*和超平面位移项的解b*
Figure BDA0003049521970000115
(5.6)将超平面法向量的解ω*和超平面位移项的解b*代入支撑向量机的决策函数J(f),得到训练好的支撑向量机,即支撑向量机检测器:
Figure BDA0003049521970000121
其中,sgn(·)为符号函数,f为输入的特征向量,若f是目标单元回波的特征向量,则J(f)=1,若f是非目标单元回波的特征向量,则J(f)=-1。
步骤6,利用支撑向量机检测器对悬停无人机进行特征检测。
(6.1)分别使用测试数据集中每一个慢时间域信号的特征值,组成测试特征向量f1,l
f1,l=(H1,Ef,1,mtfr,1,VARtfr,1,Eq,1),
其中,l∈[1,M1],M1为测试数据集中慢时间域回波的数量;
(6.2)将每一个测试特征向量f1,l分别输入到支撑向量机检测器,得到测试数据所属类别的判决值J(f1,l):
Figure BDA0003049521970000122
(6.3)根据J(f1,l),判断慢时间域信号所在的距离单元是否有目标存在,若J(f1,l)的值为1,则f1,l是目标单元回波的特征向量,即慢时间域信号所在的距离单元有目标存在,若J(f1,l)的值为-1,则f1,l是非目标单元回波的特征向量,即慢时间域信号所在的距离单元无目标存在。
下面结合仿真实验对本发明做进一步的描述:
1.仿真实验条件:
仿真实验的硬件平台为:处理器为Intel(R)Core(TM)i5-6500 CPU,主频为3.2GHz,内存8GB。
仿真实验的软件平台为:Windows 7操作系统和MATLAB R2019a。
仿真实验的雷达参数为:雷达发射信号的波长λ=0.03124m,雷达接收机的采样频率fs=100MHz,单帧数据的相参处理时间为CPI=100ms。
仿真实验的目标参数为:目标为一悬停的四旋翼无人机,无人机到雷达的初始距离R0=500m,目标速度v=0m/s,方位角α=0°,俯仰角β=0°,无人机旋翼转速fΩ=100r/s,无人机单个旋翼长度L=0.12m,无人机金属转子半径r=0.01m,无人机旋翼和金属转子的初始旋转角[θ0123]分别为区间[0,2π]上的四个随机数。
仿真实验的杂波参数为:杂波的分布为韦布尔分布,韦布尔分布的形状参数为p=5,韦布尔分布的尺度参数为q=450,杂波的功率谱为立方谱,立方谱的转折频率为fc=9Hz。
2.仿真内容及结果分析:
仿真1,在雷达脉冲重复频率为2kHz,信噪比为10dB,信杂比为0dB的条件下,仿真产生低重频情况下的悬停无人机回波,结果图2所示,其中:
图2(a)为信号时域波形图,
图2(b)为信号频谱图,
图2(c)为信号时频谱图,
图2(d)为信号倒谱图。
由图2可以看出,在脉冲重复频率较低时,悬停无人机回波中仍有明显的微多普勒特征,并且图中的微多普勒特征是由金属转子产生的。因此,可以证明根据无人机金属转子的微多普勒特征能实现低重频情况下的悬停无人机特征检测。
仿真2,在信噪比为10dB,雷达脉冲重复频率分别为20kHz和2kHz的条件下,在区间[-40dB,20dB]上,每隔5dB取一个值作为信杂比,在每一个信杂比下仿真产生1000个目标单元回波和1000个非目标单元回波,用本发明的方法分别在每一个脉冲重复频率和信杂比的组合下对仿真产生的1000个目标单元回波和1000个非目标单元回波进行特征检测,检测结果如图3所示。
从图3可以看出,当雷达脉冲重复频率为20kHz时,本发明在信杂比大于等于-25dB时能取得大于90%的检测概率,在信杂比大于等于-25dB时能取得小于10%的虚警概率;当雷达脉冲重复频率为2kHz时,本发明在信杂比大于等于-20dB时能取得大于90%的检测概率,在信杂比大于等于-25dB时能取得小于10%的虚警概率。
仿真实验结果证明,本发明在低重频、强杂波环境下能对悬停无人机进行有效的检测。

Claims (9)

1.一种基于微多普勒效应的悬停无人机特征检测方法,其特征在于,包括如下:
(1)对悬停无人机的雷达回波信号在快时间域进行下变频和脉冲压缩处理,并提取每一个距离单元的慢时间域信号;
(2)将慢时间域数据分为目标单元回波和非目标单元回波两部分,取一半目标单元回波和一半非目标单元回波组成训练数据集,将剩余的回波数据作为测试数据集;
(3)根据训练数据集和测试数据集的慢时间域数据,在不同变换域中计算表征无人机金属转子和旋翼微多普勒特性的特征值:
(3a)在时域计算训练数据的时域赫斯特指数H0和测试数据的时域赫斯特指数H1
(3b)在频域计算训练数据的频谱熵值Ef,0和测试数据的频谱熵值Ef,1
(3c)在时频域计算训练数据时频脊的均值mtfr,0和方差VARtfr,0,同时计算测试数据时频脊的均值mtfr,1和方差VARtfr,1
(3d)在倒频域计算训练数据的倒谱熵值Eq,0和测试数据的倒谱熵值Eq,1
(4)用上述计算出训练数据集中的所有特征值组成训练特征向量,对支撑向量机进行训练,生成支撑向量机检测器;
(5)使用上述计算出测试数据集的所有特征值组成特征向量,输入到支撑向量机检测器,判断是否有目标存在,实现对悬停无人机的特征检测。
2.根据权利要求1所述的方法,其特征在于,(1)中对悬停无人机的雷达回波信号在快时间域依次进行下变频和脉冲压缩处理,实现如下:
(1a)利用下变频公式,在快时间域内下变频雷达接收的悬停无人机回波,得到下变频信号h:
h=srexp(-j2πztr),
其中,sr表示雷达接收的回波信号,exp表示以自然常数e为底的指数操作,j表示虚数单位符号,z表示根据工程需求设定的雷达机动目标回波信号的载频,tr表示雷达接收回波信号快时间域的采样时间点,tr的取值范围为
Figure FDA0004076022400000021
fs表示根据工程需求设定的雷达接收机动回波信号的采样频率;
(1b)利用脉冲压缩公式,在快时间域内对下变频信号在距离维进行脉冲压缩,得到脉冲压缩信号g:
Figure FDA0004076022400000022
其中,h表示下变频信号,
Figure FDA0004076022400000023
表示卷积操作,rect(·)表示矩形窗因子,当
Figure FDA0004076022400000024
时矩形窗因子取值为0,当
Figure FDA0004076022400000025
时矩形窗因子取值为1,Tp表示雷达接收的回波信号的脉冲宽度,|·|表示取绝对值操作,k表示雷达接收的回波信号的调频斜率。
3.根据权利要求1所述的方法,其特征在于,(1)中提取出的慢时间域信号,包括非目标单元和目标单元的慢时间域回波信号,分别表示如下:
非目标单元的慢时间域回波信号x(n)为:x(n)=sc+sn,其中,n∈[1,N],N为单帧数据的脉冲个数,sc为杂波,sn为噪声;
目标单元的慢时间域回波信号x(n)为:
Figure FDA0004076022400000026
其中,nr为旋转部件的数量,a和b为常数,t=nTr为慢时间域采样点,Tr为脉冲重复时间,s0为机身回波,sR,k(t)为第k个旋翼的回波,sm,k(t)为第k个金属转子的回波。
所述的无人机旋翼回波的公式如下:
Figure FDA0004076022400000027
其中,L为旋翼的长度,exp表示以自然常数e为底的指数操作,j表示虚数单位符号,fd为无人机径向运动的多普勒频率,λ为雷达信号的波长,R0为目标无人机到雷达的初始距离,nb为单个旋翼的叶片数量,β为无人机相对雷达视线方向的俯仰角,ψt,k为物体表面的散射点P的旋转角,即雷达与旋翼所在的直线和旋翼叶片所在直线之间的夹角。
所述的无人机金属转子的回波如下:
Figure FDA0004076022400000031
其中,d为常数,J0(·)为第一类0阶贝塞尔函数,r为金属转子的半径,γt,k为金属转子表面的某散射点P所在切线和雷达与旋翼所在直线的夹角,ψt,k为物体表面的散射点P的旋转角,U1(t)和U2(t)为周期矩形脉冲信号:
Figure FDA0004076022400000032
U2(t)=1-U1(t)。
4.根据权利要求1所述的方法,其特征在于,(3a)中在时域计算训练数据的时域赫斯特指数H0和测试数据的时域赫斯特指数H1,实现如下:
(3a1)分别从训练数据集中取慢时间域信号x0(n),计算训练数据的累积离差Y0(n),从测试数据集中取慢时间域信号x1(n),计算测试数据的累积离差Y1(n):
Figure FDA0004076022400000033
Figure FDA0004076022400000034
其中,n∈[1,N],N为单帧数据的脉冲个数,<·>表示取均值的操作;
(3a2)分别从训练数据的累积离差序列Y0(n)的开头和末尾开始划分出
Figure FDA0004076022400000041
个连续且长度为m的短序列,构成2Nm个短序列,同时,对测试数据的累积离差序列Y1(n)执行同样的划分操作,也构成2Nm个短序列;
(3a3)使用最小二乘法分别对序列Y0(n)和序列Y1(n)划分出的短序列进行拟合,分别得到训练数据的拟合多项式y0(i)和测试数据的拟合多项式y1(i),其中,i的取值范围为[1,2Nm];
(3a4)分别对训练数据和测试数据进行消除趋势估计,获得训练数据的拟合方差Yυ,0(s,m)和测试数据的拟合方差Yυ,1(s,m):
Figure FDA0004076022400000042
Figure FDA0004076022400000043
(3a5)对(3a4)得到的两个拟合方差分别取均值,再开方,得到训练数据的二阶波动函数F0(m)和测试数据的二阶波动函数F1(m):
Figure FDA0004076022400000044
Figure FDA0004076022400000045
其中,s为拟合多项式的序号;
(3a6)根据训练数据的二阶波动函数F0(m)和测试数据的二阶波动函数F1(m),使用最小二乘法,分别计算训练数据的时域赫斯特指数H0和测试数据的时域赫斯特指数H1
Figure FDA0004076022400000046
Figure FDA0004076022400000051
5.根据权利要求1所述的方法,其特征在于,(3b)中在频域计算训练数据的频谱熵值Ef,0和测试数据的频谱熵值Ef,1,实现如下:
(3b1)分别对训练数据的频谱X0(f)和测试数据的频谱X1(f)进行归一化,得到训练数据的归一化频谱X0,1(f)和测试数据的归一化频谱X1,1(f):
Figure FDA0004076022400000052
Figure FDA0004076022400000053
其中,fr为雷达的脉冲重复频率;
(3b2)分别计算训练数据的频谱熵值Ef,0和测试数据的频谱熵值Ef,1
Figure FDA0004076022400000054
Figure FDA0004076022400000055
6.根据权利要求1所述的方法,其特征在于,(3c)中在时频域计算训练数据时频脊的均值mtfr,0和方差VARtfr,0,同时计算测试数据时频脊的均值mtfr,1和方差VARtfr,1,实现如下:
(3c1)将短时傅里叶变换的短时窗设置为汉明窗,为使金属转子的微多普勒特征易于观察和提取,设窗长应大于旋翼旋转周期的2倍,即:
Lw>2T
其中,Lw为短时傅里叶变换的窗长,T为无人机旋翼和金属转子的旋转周期;
(3c2)分别对训练数据x0(n)和测试数据x1(n)作短时傅里叶变换,得到训练数据的时频谱STFTx0(n,f)和测试数据的时频谱STFTx1(n,f):
Figure FDA0004076022400000061
Figure FDA0004076022400000062
其中,n∈[1,N],N为单帧数据的脉冲个数,w(·)为短时窗,exp表示以自然常数e为底的指数操作,j表示虚数单位符号;
(3c3)分别对训练数据和测试数据的时频谱进行截取,得到去除机身信号的训练数据时频谱STFTx0,0(n,f)和去除机身信号的测试数据时频谱STFTx1,0(n,f):
Figure FDA0004076022400000063
Figure FDA0004076022400000064
(3c4)根据去除机身信号的时频谱,分别计算为训练数据中无人机金属转子的时频脊R0(n)和测试数据中无人机金属转子的时频脊R1(n):
R0(n)={|fm||STFTx0,0(n,fm)=max[STFTx0,0(n,f)]}
R1(n)={|fm||STFTx1,0(n,fm)=max[STFTx1,0(n,f)]}
(3c5)根据(3c4)的结果,分别计算训练数据时频脊的均值mtfr,0和测试数据时频脊的均值mtfr,1
Figure FDA0004076022400000065
Figure FDA0004076022400000066
(3c6)根据(3c4)和(3c5)的结果,分别计算训练数据时频脊的方差VARtfr,0和测试数据时频脊的方差VARtfr,1
Figure FDA0004076022400000071
Figure FDA0004076022400000072
7.根据权利要求1所述的方法,其特征在于,(3d)中在倒频域计算训练数据的倒谱熵值Eq,0和测试数据的倒谱熵值Eq,1,实现如下:
(3d1)分别对训练数据和测试数据的倒谱进行截取,得到去除机身信号的训练数据倒谱C0,0(q)和去除机身信号的测试数据倒谱C1,0(q):
C0,0(q)=C0(q),q∈(T/3,NTr/2)
C1,0(q)=C1(q),q∈(T/3,NTr/2)
其中,C0(q)为训练数据的倒谱,C1(q)为测试数据的倒谱,T为旋翼和金属转子的旋转周期,N为单帧数据的脉冲数,Tr为脉冲重复周期;
(3d2)对截取的倒谱进行归一化,分别计算训练数据的归一化倒谱C0,1(q)和测试数据的归一化倒谱C1,1(q):
Figure FDA0004076022400000073
Figure FDA0004076022400000074
(3d3)根据(3d2)的结果,分别计算训练数据的倒谱熵值Eq,0和测试数据的倒谱熵值Eq,1
Figure FDA0004076022400000075
Figure FDA0004076022400000076
8.根据权利要求1所述的方法,其特征在于,(4)中利用训练特征向量,对支撑向量机进行训练,实现如下:
(4a)为训练特征向量标注类别,得到训练样本集Ω0
Ω0={(f0,i,yi)|i=1,2,...,M}
其中,M为训练数据集中慢时间域回波的数据量,f0,i为训练特征向量,yi为回波信号的类别,若回波来自目标单元,则yi=1,若回波来自非目标单元,则yi=-1;
(4b)根据训练样本集构造一个最优超平面作为决策曲面,使得目标单元回波和非目标单元回波样本之间的隔离边缘最大化,最优超平面g(f)满足如下条件:
yig(f)=yiTf0,i+b)≥1,i=1,2,...,M
其中,ω为超平面的法向量,b为超平面的位移项,(·)T为转置操作;
(4c)根据(4b)中的限制条件,将最优超平面参数的求解问题转化为如下的二次规划问题:
Figure FDA0004076022400000081
再通过拉格朗日乘数法,将上述二次规划转化为对偶问题Q(α):
Figure FDA0004076022400000082
其中,α=(α1;...;αi;...;αM)为拉格朗日系数,αi为拉格朗日系数α中的第i个值,αi≥0,i∈[1,M];
(4d)对(4c)中的对偶问题Q(α)进行求解,得到拉格朗日系数α的解为
Figure FDA0004076022400000083
再由拉格朗日系数的解α*计算超平面法向量的解ω*和超平面位移项的解b*
Figure FDA0004076022400000091
(4e)将超平面法向量的解ω*和超平面位移项的解b*代入支撑向量机的决策函数J(f),得到训练好的支撑向量机,即支撑向量机检测器:
Figure FDA0004076022400000092
其中,sgn(·)为符号函数,f为输入的特征向量,若f是目标单元回波的特征向量,则J(f)=1,若f是非目标单元回波的特征向量,则J(f)=-1;
9.根据权利要求1所述的方法,其特征在于,(5)中将测试数据集的特征向量,输入到支撑向量机检测器,判断是否有目标存在,是使用支撑向量机检测器对测试数据的特征向量f1进行计算,得到计算结果J(f1),若J(f1)=1,则f1是目标单元回波的特征向量,即测试数据对应的距离单元有目标存在,若J(f1)=-1,则f1是非目标单元回波的特征向量,即测试数据对应的距离单元无目标存在。
CN202110481764.7A 2021-04-30 2021-04-30 基于微多普勒效应的悬停无人机特征检测方法 Active CN113238210B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110481764.7A CN113238210B (zh) 2021-04-30 2021-04-30 基于微多普勒效应的悬停无人机特征检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110481764.7A CN113238210B (zh) 2021-04-30 2021-04-30 基于微多普勒效应的悬停无人机特征检测方法

Publications (2)

Publication Number Publication Date
CN113238210A CN113238210A (zh) 2021-08-10
CN113238210B true CN113238210B (zh) 2023-04-07

Family

ID=77131771

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110481764.7A Active CN113238210B (zh) 2021-04-30 2021-04-30 基于微多普勒效应的悬停无人机特征检测方法

Country Status (1)

Country Link
CN (1) CN113238210B (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111398909A (zh) * 2020-03-09 2020-07-10 哈尔滨工业大学 一种基于倒谱分析的杂波环境无人机检测方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109633629B (zh) * 2018-10-26 2020-11-03 上海无线电设备研究所 太赫兹频段单旋翼无人机目标特性微多普勒特征提取方法
EP3715901A1 (en) * 2019-03-29 2020-09-30 Robin Radar Facilities BV Detection and classification of unmanned aerial vehicles
US11442159B2 (en) * 2019-07-19 2022-09-13 The Regents Of The University Of California Multi-spectral THz micro-doppler radar based on silicon-based picosecond pulse radiators

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111398909A (zh) * 2020-03-09 2020-07-10 哈尔滨工业大学 一种基于倒谱分析的杂波环境无人机检测方法

Also Published As

Publication number Publication date
CN113238210A (zh) 2021-08-10

Similar Documents

Publication Publication Date Title
CN111352102B (zh) 一种基于调频连续波雷达的多目标个数检测方法及装置
CN103197301B (zh) 海面微动目标Radon-线性正则变换长时间相参积累检测方法
CN103176178B (zh) 雷达动目标Radon-分数阶傅里叶变换长时间相参积累检测方法
CN102788969B (zh) 基于短时分数阶傅里叶变换的海面微动目标检测和特征提取方法
CN109521410B (zh) 基于时间反转变换的高速机动目标相参积累检测方法
CN107450055B (zh) 基于离散线性调频傅立叶变换的高速机动目标检测方法
CN112162273B (zh) 一种基于奇异矢量的多旋翼无人机物理参数提取方法
CN111736128A (zh) 基于skt-siaf-mscft的相参积累方法
CN108196241B (zh) 一种基于Hough变换的高速动目标速度估计方法
CN112198487B (zh) 一种风电场杂波背景下的目标检测方法
Molchanov et al. On micro-Doppler period estimation
Cosoli et al. A real-time and offline quality control methodology for SeaSonde high-frequency radar currents
CN113238210B (zh) 基于微多普勒效应的悬停无人机特征检测方法
Ivanov et al. CFAR multi-target detection based on non-central Chi-square distribution for FMCW
CN114942419B (zh) 长积累时间下的船只散射点三自由度微动特征提取方法
CN115932824A (zh) 一种基于多天线的fmcw雷达测距方法及系统
CN113009443B (zh) 一种基于图连通密度的海面目标检测方法及其装置
CN112906476B (zh) 一种基于信杂噪比损失的机载雷达训练样本选择方法
CN114740445A (zh) 一种基于平均谱半径的海面漂浮小目标检测方法
CN113820703A (zh) 一种基于散射变换的无人机目标旋翼参数估计方法
Bai et al. Floating small target detection based on the dual-polarization cross-time-frequency distribution in sea clutter
Zhang et al. Target classification with low-resolution radars based on multifractal features in fractional fourier domain
Xie et al. Micro-Doppler-based method for rotor parameter estimation under short dwell time conditions
CN114578308B (zh) 一种基于混合多普勒的旋翼目标特征提取方法
CN113702969B (zh) 基于自适应stft方法的微多普勒信号参数估计方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant