基于GLCT-GPTF的多旋翼目标参数估计方法

基于glct
‑
gptf的多旋翼目标参数估计方法
技术领域
1.本发明属于雷达技术领域，更进一步涉及雷达目标识别技术领域中的一种基于广义线性调频小波变换
‑
广义参数化时频分析glct
‑
gptf(general linear chirplet transform
‑
general parameterized time
‑
frequency)的多旋翼目标参数估计方法。本发明可用于对直升飞机和螺旋桨飞机等具有多个旋翼的目标进行旋翼物理参数估计，所估计的旋翼物理参数可用于多旋翼目标的识别。


背景技术：

2.旋翼桨叶的旋转会对雷达的电磁回波产生调制作用，在回波的时频谱图中表现为正弦形式的微多普勒曲线，微多普勒曲线的周期与旋翼的转动周期一致，最大幅度由雷达波长、旋翼桨叶长度和旋翼转速决定。由于不同类型的多旋翼目标的物理参数不同，因此可以从多旋翼目标的雷达回波中估计旋翼物理参数，将估计得到的旋翼物理参数与旋翼目标知识库中旋翼物理参数相比较，可以实现多旋翼目标的分类。现有的方法主要可以分为两类。第一类方法是非参数化方法，这类方法主要是在时频谱图的基础上通过图像形态学处理、峰值检测等方法提取桨叶的微多普勒展宽，结合微动周期估计方法得到旋翼转速，最后根据旋翼桨叶长度与微多普勒展宽和旋翼转动频率之间的关系实现旋翼桨叶长度的估计。上述非参数方法的性能均受到时频分析方法分辨率的限制，旋翼物理参数估计精度不高。参数化方法将旋翼物理参数估计问题转化为时频谱图中的正弦微多普勒曲线提取问题，如霍夫变换、逆拉东变换等，然而由于多旋翼目标时频图中多条微多普勒曲线的存在，回波在时频谱图上的能量集中度较差，降低了旋翼物理参数估计精度。
3.电子科技大学在其申请的专利文献“一种基于奇异矢量的多旋翼无人机物理参数提取方法”(申请号202011024462.9，申请公布号cn 112162273 a)中公开了一种基于奇异矢量的多旋翼无人机物理参数提取方法。该方法的具体步骤如下：对多旋翼无人机雷达回波信号进行小波分解，提取叶片旋转部件产生的低频分量，然后，将此分量进行短时傅里叶变化得到时频谱图，由时频谱图奇异值分解后最大奇异值对应的奇异矢量提取旋翼的转速和桨叶长度等物理参数。该方法存在的不足之处是：使用上述方法对参数进行估计时，基于短时傅里叶变换的时频谱图能量集中度不高，降低了奇异矢量提取的准确性以及旋翼物理参数估计的精度。
4.电子科技大学在其申请的专利文献“一种直升机旋翼物理参数提取方法”(申请号201910519253.2，申请公布号cn 110133600 a)中公开了一种基于时频谱图的直升机旋翼物理参数提取方法。该方法的具体步骤如下：首先对直升机的时频谱图进行滤波，将滤波后的时频谱图按照经验得到的两个阈值进行三值量化和分割，减弱时频谱图的背景噪声，提升图像清晰度，以准确提取时频信号线；再利用最小二乘法估计出直升机旋翼的旋转周期；最后通过频谱宽度与桨叶长度的关系，估算出桨叶的长度。该方法存在的不足之处是：使用上述方法对参数进行估计时时频谱图量化时的门限是根据经验得到的，对于不同类型的多旋翼目标，此经验门限很难预先确定，实际中对于不同类型的多旋翼目标难以直接应用。


技术实现要素：

5.本发明的目的在于针对上述现有技术的不足，提出一种基于glct
‑
gptf多旋翼目标物理参数估计方法，旨在解决现有技术的基于短时傅里叶变换时频谱图的旋翼物理参数估计方法时频谱图能量集中度不高以及需要预先设置时频谱图量化门限影响旋翼物理参数估计精度的问题。
6.实现本发明目的的思路是：本发明通过由变换核参数向量决定的glct
‑
gptf变换旋转平移多旋翼目标雷达回波信号，利用天牛群优化算法获得与多旋翼目标雷达回波信号匹配的变换核参数向量，使得旋转平移后信号趋于平稳，能够有效提升时频谱图能量集中度。由于本发明基于glct
‑
gptf变换的信号平稳度衡量变换核参数向量与多旋翼目标雷达回波信号的匹配度，将对旋翼参数估计的问题，转化为旋转平移后信号的平稳度最大化问题，变换核参数向量中旋翼的最大多普勒频率是由旋翼的转速和旋翼的桨叶长度共同决定的，通过变换核参数向量与旋翼物理参数之间的映射关系可以直接实现旋翼转速以及桨叶长度的估计，由此无需通过时频谱图量化得到最大多普勒频率的估计，避免了量化门限的预先设置而影响旋翼物理参数估计精度的问题。
7.实现本发明目的的具体步骤如下：
8.步骤1，计算每个旋翼目标的转速：
9.(1a)计算雷达回波信号在每个时刻和每个延迟时间的三阶瞬时矩；
10.(1b)计算雷达回波信号在每个频率和每个延迟时间的三阶模糊值；
11.(1c)按照下式，计算雷达回波信号在每个延迟时间的三阶模糊熵：
[0012][0013]
其中，ε(τ)表示雷达回波信号在延迟τ时间的三阶模糊熵，0＜τ≤n，n表示雷达回波信号中的时刻总数,n≥2000，||表示取模操作，c(f,τ)表示雷达回波信号在第f个频率，延迟τ时间的三阶模糊值，1≤f≤n，ln()表示自然常数e为底的对数操作；
[0014]
(1d)按照下式，计算旋翼的转速：
[0015][0016]
其中，ω表示旋翼的转速，π表示圆周率，prf表示雷达回波信号脉冲重复频率，prf≥10000hz,n
l
表示旋翼桨叶的总数，s表示符号系数，当n
l
为偶数时，s＝1，当n
l
为奇数时，s＝2，t表示在所有延迟时间的模糊熵的最大值对应的延迟时间；
[0017]
步骤2，利用广义线性调频小波变换
‑
广义参数化时频分析glct
‑
gptf变换计算旋转平移后信号的平稳度：
[0018]
(2a)按照下式，构建与多旋翼雷达回波信号相匹配的正弦形式的变换核：
[0019][0020]
其中，κ
p
(t)表示第t个时刻，基于由旋翼的最大多普勒频率、旋翼的转速、旋翼的初始相位构成的变换核参数向量p得到的正弦形式的变换核，p(1)、p(2)、p(3)分别表示变换核参数向量p中的第一个元素，第二个元素和第三个元素；
[0021]
(2b)按照下式，计算三维时频谱图中每个像素的模值：
[0022][0023]
其中，s(t,f,c,p)表示第t个时刻，第f个频率，第c个调频率下，基于变换核参数向量p得到的三维时频谱图的模值，1≤c≤5，s(τ)表示雷达回波信号向量中第τ个时刻的雷达回波信号，表示第τ个时刻由变换核参数向量p得到的信号旋转算子，exp表示以自然常数e为底的指数操作，j表示虚数符号，κ
p
(τ)表示第τ个时刻由变换核参数向量p得到正弦形式的变换核，sin表示正弦操作，表示第τ个时刻由变换核参数向量p得到的信号平移算子，κ
p
(t)表示第t个时刻由变换核参数向量p得到的正弦形式的变换核，tan表示正切操作；
[0024]
(2c)按照下式，计算旋转平移后信号的平稳度：
[0025][0026]
其中，f表示旋转平移后信号的平稳度，e表示由变换核所拟合的曲线在二维时频谱图s
c
(t,f,p)上对应像素的模值的总和，(t,f,p)上对应像素的模值的总和，表示三维时频谱图上第t个时刻，第f个频率下的最优调频率，c表示由变换核所拟合的曲线在三维时频谱图上对应像素处最优调频率不等于3的点数的总和；
[0027]
步骤3，利用天牛群优化算法获得与多旋翼目标雷达回波信号匹配的变换核参数向量：
[0028]
使用天牛群优化算法迭代更新变换核参数向量p中旋翼的最大多普勒频率和旋翼的初始相位，直到信号的平稳度不再变化为止，将迭代终止时天牛群中所有天牛的信号平稳度最大对应的变换核参数向量，作为与多旋翼目标雷达回波信号相匹配的变换核参数向量；
[0029]
步骤4，计算旋翼的桨叶长度。
[0030]
本发明与现有技术相比，具有以下优点：
[0031]
第一，本发明利用广义线性调频小波变换
‑
广义参数化时频分析glct
‑
gptf变换计算旋转平移后信号的平稳度，克服了现有技术中基于短时傅里叶变换时频谱图的旋翼物理参数估计方法时频谱图能量集中度不高的缺陷，使得本发明具有更精确的刻画多旋翼目标雷达回波信号时频特性的优点。
[0032]
第二，本发明利用天牛群优化算法获得与多旋翼目标雷达回波信号匹配的变换核参数向量，克服了现有技术需要预先设置时频谱图量化门限影响旋翼物理参数估计精度的问题，使得本发明提高了旋翼的转速以及桨叶长度的估计精度。
附图说明
[0033]
图1为本发明的流程图。
具体实施方式
[0034]
下面结合附图1，对本发明的具体步骤做进一步的描述。
[0035]
步骤1，计算每个旋翼目标的转速。
[0036]
按照下式，计算雷达回波信号在每个时刻和每个延迟时间的三阶瞬时矩：
[0037]
c(t,τ)＝s(t)[s(t+τ)]
*
[s(t+τ)]
*
s(t+2τ)
[0038]
其中，c(t,τ)表示雷达回波信号中第t个时刻，延迟τ时间的三阶瞬时矩，1≤t≤n，s(t)表示雷达回波信号向量中第t个时刻的雷达回波信号，s表示雷达回波信号向量，s(t+τ)表示雷达回波信号向量中第t+τ个时刻的雷达回波信号，*表示共轭操作，s(t+2τ)表示雷达回波信号向量中第t+2τ个时刻的雷达回波信号。
[0039]
按照下式，计算雷达回波信号在每个频率和每个延迟时间的三阶模糊值：
[0040]
c(f,τ)＝∫c(t,τ)exp(
‑
j2πft)dt
[0041]
其中，c(f,τ)表示雷达回波信号在第f个频率，延迟τ时间的三阶模糊值。
[0042]
按照下式，计算雷达回波信号在每个延迟时间的三阶模糊熵：
[0043][0044]
其中，ε(τ)表示雷达回波信号在延迟τ时间的三阶模糊熵，0＜τ≤n，n表示雷达回波信号中的时刻总数,n≥2000，||表示取模操作，c(f,τ)表示雷达回波信号在第f个频率，延迟τ时间的三阶模糊值，1≤f≤n，ln()表示自然常数e为底的对数操作。
[0045]
按照下式，计算旋翼的转速：
[0046][0047]
其中，ω表示旋翼的转速，π表示圆周率，prf表示雷达回波信号脉冲重复频率，prf≥10000hz,n
l
表示旋翼桨叶的总数，s表示符号系数，当n
l
为偶数时，s＝1，当n
l
为奇数时，s＝2，t表示在所有延迟时间的模糊熵的最大值对应的延迟时间。
[0048]
步骤2，利用广义线性调频小波变换
‑
广义参数化时频分析glct
‑
gptf变换计算旋转平移后信号的平稳度。
[0049]
按照下式，构建与多旋翼雷达回波信号相匹配的正弦形式的变换核：
[0050][0051]
其中，κ
p
(t)表示第t个时刻，基于由旋翼的最大多普勒频率、旋翼的转速、旋翼的初始相位构成的变换核参数向量p得到的正弦形式的变换核，p(1)、p(2)、p(3)分别表示变换核参数向量p中的第一个元素，第二个元素和第三个元素。
[0052]
按照下式，计算三维时频谱图中每个像素的模值：
[0053][0054]
其中，s(t,f,c,p)表示第t个时刻，第f个频率，第c个调频率下，基于变换核参数向量p得到的三维时频谱图的模值，1≤c≤5，s(τ)表示雷达回波信号向量中第τ个时刻的雷达回波信号，表示第τ个时刻由变换核参数向量p得到的信号旋转算子，exp表示以自然常数e为底的指数操作，j表示虚数符号，κ
p
(τ)表示第τ个时刻由变换核参数向量p得到正弦形式的变换核，sin表示正弦操作，表示第τ个时刻由变换核参数向量p得到的信号平移算子，κ
p
(t)表示第t个时刻由变换核参数向量p得到的正弦形式的变换核，tan表示正切操作。
[0055]
所述的三维时频谱图上第t个时刻，第f个频率下的最优调频率是由下式得到的：
[0056][0057]
其中，表示给定变换核参数p时第t个时刻，第f个频率下的三维时频谱图的模值s(t,f,c,p)取得最大值时c的取值。
[0058]
按照下式，计算旋转平移后信号的平稳度：
[0059][0060]
其中，f表示旋转平移后信号的平稳度，e表示由变换核所拟合的曲线在二维时频谱图s
c
(t,f,p)上对应像素的模值的总和，(t,f,p)上对应像素的模值的总和，表示三维时频谱图上第t个时刻，第f个频率下的最优调频率，c表示由变换核所拟合的曲线在三维时频谱图上对应像素处最优调频率不等于3的点数的总和。
[0061]
步骤3，利用天牛群优化算法获得与多旋翼目标雷达回波信号匹配的变换核参数向量。
[0062]
使用天牛群优化算法迭代更新变换核参数向量p中旋翼的最大多普勒频率和旋翼的初始相位，直到信号的平稳度不再变化为止，将迭代终止时天牛群中所有天牛的信号平稳度最大对应的变换核参数向量，作为与多旋翼目标雷达回波信号相匹配的变换核参数向量。
[0063]
天牛群优化算法是粒子群优化算法和天牛须算法的结合，可以比较好地避免传统粒子群算法存在的稳定性不高、运算量较大、容易陷入局部最优的问题。粒子群算法中的每个粒子都被替换为一只天牛同时进行搜索，天牛个体的位置和速度的初始化由下式得到：
[0064][0065]
[0066]
其中，x
im
表示第i次迭代时第m个天牛的位置向量，0＜i≤300，0＜m≤30，表示与旋翼的最大多普勒频率对应的第i次迭代时第m个天牛的位置向量的第一个元素，从区间中随机选取，表示第i次迭代时第m个天牛的位置向量的第二个元素，从区间[0,2π]中随机选取，v
im
表示第i次迭代时第m个天牛的速度向量，表示第i次迭代时第m个天牛的第一维的速度，表示第i次迭代时第m个天牛的第二维的速度，和均从区间[
‑
2,2]中随机选取。
[0067]
在迭代过程中，每个天牛的位置的更新规则不但结合所有天牛的最优全局位置和该天牛个体的自身最优位置，并且引入了天牛触须的搜索机制，在迭代过程中对天牛个体左、右两侧触须位置向量对应的信号平稳度进行比较，将平稳度更大的位置用于更新天牛个体的位置。利用天牛群优化算法可以有效提升变换核参数向量估计的精度以及速度。
[0068]
步骤4，按照下式，计算旋翼的桨叶长度：
[0069][0070]
其中，l表示旋翼的桨叶长度，λ表示雷达回波信号的波长，p
max
表示与多旋翼目标雷达回波信号匹配的变换核参数向量，p
max
(1)表示与多旋翼目标雷达回波信号匹配的变换核参数向量中的第一个元素。
[0071]
下面结合仿真实验对本发明的效果做进一步的描述。
[0072]
1.仿真实验条件：
[0073]
本发明的仿真实验的硬件平台为：处理器为intel i7 6700k cpu，主频为4.0ghz，内存16gb。
[0074]
本发明的仿真实验的软件平台为：windows 10操作系统和matlab2020a。
[0075]
本发明仿真实验所使用的雷达回波数据是在合作情况下采集到的“运八”飞机的高信噪比实测雷达数据。“运八”飞机为四旋翼螺旋桨飞机，其桨叶长度l＝1.8m，旋翼转动速度ω
r
＝40πrad/s，每个旋翼的桨叶数目n＝4；雷达脉冲重复频率12khz，驻留时间为1s，载频为5ghz。
[0076]
2.仿真内容及其结果分析：
[0077]
本发明仿真实验是采用本发明的方法和两个现有技术(rsp
‑
cfd方法、瞬时频率
‑
fft方法)分别对“运八”飞机的高信噪比实测雷达数据的旋翼转速和桨叶长度进行估计。
[0078]
在仿真实验中，采用的两个现有技术是指：
[0079]
现有技术rsp
‑
cfd方法是指，何炜琨等人在“基于rsp
‑
cfd方法的小型旋翼无人机微动特征提取[j].信号处理,2021,37(3):399
‑
408”中提出的一种多旋翼无人机目标旋翼参数估计方法，简称rsp
‑
cfd方法。
[0080]
现有技术瞬时频率
‑
fft方法是指，马娇等人在“多旋翼无人机微多普勒特性分析与特征提取[j].中国科学院大学学报,2019,36(2):235
‑
243”中提出的一种多旋翼无人机目标旋翼参数估计方法，简称瞬时频率
‑
fft方法。
[0081]
对本发明的仿真结果得到的旋翼的最大多普勒频率、旋翼的转速以及旋翼的桨叶长度评价结果如表1所示。表1中本发明“最大多普勒频率”的参数是本发明步骤3中通过天
牛群优化算法获得的与多旋翼目标雷达回波信号匹配的变换核参数向量中的第一个元素。表1中本发明“转速”的参数是利用本发明步骤1中三阶模糊熵估计得到的旋翼的转速。表1中本发明“桨叶长度”的参数是利用本发明步骤4中旋翼的桨叶长度的计算公式得到。表1中rsp
‑
cfd方法“最大多普勒频率”的参数是在雷达回波信号的量化后的rsp时频谱图上通过极大值搜索得到。表1中rsp
‑
cfd方法“转速”的参数是在雷达回波信号的量化后的cfd时频谱图上通过极大值搜索得到。表1中rsp
‑
cfd方法“桨叶长度”的参数是利用最大多普勒频率与转速以及桨叶长度之间的映射关系计算得到。表1中瞬时频率
‑
fft方法“最大多普勒频率”的参数是在雷达回波信号的量化后的短时傅里叶变换时频谱图上通过极大值搜索得到。表1中瞬时频率
‑
fft方法“转速”的参数是对雷达回波信号的量化后的短时傅里叶变换时频谱图进行傅里叶变换搜索极大值得到。表1中瞬时频率
‑
fft方法“桨叶长度”的参数是利用最大多普勒频率与转速以及桨叶长度之间的映射关系计算得到。表1中真实值“转速”的参数以及“桨叶长度”的参数是由合作方仪器测量得到。表1中真实值“最大多普勒频率”的参数是利用最大多普勒频率与转速以及桨叶长度之间的映射关系计算得到。
[0082]
表1 仿真实验中本发明和各现有技术旋翼参数估计结果表
[0083] 真实值本发明方法rsp
‑
cfd瞬时频率
‑
fft最大多普勒频率(hz)4484.854411.094210.204250.50转速(πrad/s)40.0040.0039.9239.96桨叶长度(m)1.801.771.691.70
[0084]
由表1可以看出，本发明方法的旋翼的最大多普勒频率、转速以及桨叶长度的估计值与真实值的误差均低于2种现有技术方法，证明本发明具有更高的旋翼参数估计精度。
[0085]
以上仿真实验表明：本发明方法利用三阶模糊熵实现了旋翼转速的精确估计，通过由变换核参数向量决定的glct
‑
gptf变换旋转平移多旋翼目标雷达回波信号，利用天牛群优化算法获得与多旋翼目标雷达回波信号匹配的变换核参数向量，使得旋转平移后信号趋于平稳，能够有效提升时频谱图能量集中度，基于glct
‑
gptf变换得到的时频谱图将旋翼参数估计问题转化为旋转平移后信号的平稳度最大化问题，无需通过时频谱图量化得到最大多普勒频率的估计，解决了现有技术方法中存在的时频谱图能量集中度不高的问题以及需要预先设置时频谱图量化门限而影响旋翼物理参数估计精度的问题，是一种非常实用的多旋翼目标参数估计方法。