CN107622036B - 一种基于蚁群优化的多项式相位信号自适应时频变换方法 - Google Patents

Abstract

本发明提供一种基于蚁群优化的多项式相位信号自适应时频变换方法，能够完成多项式相位信号的时频分解，其中分解得到的每一个信号分量都为任一时刻都只对应一个频点的单分量，然后利用各信号分量、各时刻的瞬时频率取值，通过仅保留主瓣响应的Sinc函数直接计算生成相应时刻对应的信号频率分布，克服了传统的时频变换中一个时刻对应多个频点的非单分量存在交叉项的缺陷，最终输出无任何交叉项干扰且时频联合分辨率较优的时频分布；本发明原理简单，操作方便，可有效克服经典时频分析方法交叉项干扰的不利影响以及时频联合分辨率的损失，能够有效提升非平稳多项式相位信号时频分析的质量和效益。

Description

一种基于蚁群优化的多项式相位信号自适应时频变换方法

技术领域

本发明属于信号处理领域，尤其涉及一种基于蚁群优化的多项式相位信号自适应时频变换方法。

背景技术

许多天然和人工的信号，譬如语音、生物医学信号、在色散媒质中传播的波、机械振动、动物叫声，音乐、雷达、声纳信号等，都是典型的非平稳信号，其特点是持续时间有限，并且频率是时变的,具有非平稳、非线性、非均匀、非结构、非确定、非可积、非可逆、非晶态、非规则、非连续、非光滑、非周期、非对称等特点。时频联合分析(joint time-frequencyanalysis，简称时频分析)正是着眼于真实信号组成成分的时变特征，将一个一维的时间信号以二维的时间-频率密度函数形式表示出来，旨在揭示信号中包含了多少频率分量，以及每一分量是怎样随时间变化的。

1948年，法国学者J.Ville将匈牙利布达佩斯出生的美籍物理学家E.P.Wigner在1932年提出的Wigner分布引入信号处理领域，得到了称为“Wigner-Ville分布”(Wigner-Ville distribution，WVD)。后续学者起而效仿，提出了一些新型的时频分布。整个时频分析的历史，几乎就是一部与WVD的不足作斗争的历史。按照各派的本质特征，可将形形色色的时频分布归入如下几类：(1)线性时频表示；(2)Cohen类双线性时频分布；(3)仿射类双线性时频分布；(4)重排类双线性时频分布；(5)自适应核函数类时频分布；(6)参数化时频分布。

其中，线性时频变换Gabor变换、STFT的时频分辨率受制于窗函数的形状和宽度。小波分析本质上是一种时间-尺度分析，更适于分析具有自相似结构的信号(如分形)和突变(瞬态)信号，而从刻画信号的时变结构角度看，小波变换的结果往往难于解释。

Cohen类双线性时频分布的实质，是将信号的能量(信号的某种平方形式)分布于时频平面内，其基础则为WVD。但WVD不是线性的，即两信号之和的WVD并非每一个信号的WVD之和，其中多出一个附加项。交叉项对WVD影响之烈，可见一斑。交叉项是实的，混迹于自项成分之间，且幅度较大；另外，交叉项是振荡型的，每两个信号分量就会产生一个交叉项。若信号有N个分量，则会产生

个交叉项。交叉项的存在严重地干扰着人们对WVD的解释；当信号的组成成分变得复杂时，WVD给出的时频分布甚至变得毫无意义。

为了解决交叉项干扰的影响，人们相继提出了若干种不同的时频分析方法，其中Cohen类核函数时频分析方法通过设计二维核函数(二维滤波器)，产生具有所需特性的时频分布。但是，该类时频分布用平滑的方法抑制交叉项，是以牺牲整个分布的时频分辨率为代价的。

前面介绍的各种分布(小波变换除外)对信号的时间和频率局部特性的刻画，是通过时移及频移变换实现的；与此对照的是，仿射类的分布则是通过时移和伸缩变换实现的。

仿射类中最著名的分布当推尺度图(Scale gram)，即小波变换的平方。由于这一类分布的基础仍是WVD，因此，WVD自然地成为其成员之一。事实上，WVD正是连接Cohen类和仿射类的纽带。前者基于对WVD进行的时频平滑(time-frequency smoothing)，而后者则基于仿射平滑(affine smoothing)。

Cohen类和仿射类时频分布，通过对WVD进行时频平滑及时间-尺度平滑处理，如图4所示，可以极大地抑制交叉项干扰，但仍有不少交叉项残留，而且有些分布还会引入一些新的交叉项。为进一步提高这两类分布的性能，K.Kdoera等首先提出对时频平面进行重排的思想，此后F.Auger和P.Flandrin拓展并完善了重排的方法。

在Cohen类、仿射类和重排类的双线性时频分布中，每一种分布均与一固定的核函数相对应，正是该核函数决定了相应分布的交叉项抑制特性。不言而喻，一种核函数只对一类信号有效，因而所有上述三类中的双线性时频分布均缺乏对信号的适应性。

前面介绍的各种时频分析方法，均为非参数化的方法，它们均没有先验地假定信号是由何种模型信号组成的。而参数化时频分析(parametric time frequency analysis)方法，则根据对信号层次结构的分析，构造出与信号层次结构最佳匹配的信号模型，因而能浓缩信号的信息，简化信号的表示，并由此得到信号的时频分布。

在线性时频表示(原子分解)方法中，若选用的原子与信号的主要成分相似，则仅需少数原子的线性组合，就能比较精确地表示信号，分解的结果将是稀疏的(sparse)。反之，如果原子的性状与信号的主要结构相去甚远或迥异，那么就需要用大量的，甚至无穷多的原子，才能足够精确地组装成原信号，信号的信息将弥散在太多的原子上，不利于有效地表示信号。所以在采用原子分解方法时，必须根据信号的局部结构特征，自适应地选择原子的组合，以期用尽可能少的原子来分解信号。1994年S.Qian和D.Chen开创了参数化时频分析之先河(其创造性思想可追溯到1988年)，提出了“自适应展开’，(adaptive expansion)算法；1993年S.Mallat和Z.Zhang提出了自适应展开的姊妹算法—“matching pursuit”算法。这两种算法尽管名称不同，但本质上没有不可调和的异点，实质上等价。

自适应匹配投影塔形分解算法的实质，是用原子的时频能量分布逼近原信号的时频能量分布。由于S.Qian和S.Mallat采用的是频率不变的Gabor原子，因此迭代算法对时频平面的划分是一种格型分割。这种算法对时不变的频率分量效果很好，但当待分析的信号是Chirp信号时，这种匹配相当于零阶曲线逼近，势必造成分解过程存在许多截断和分量之间的混合畸变。为克服这一缺陷，S.Mann和S.Haykin等人几乎同时提出采用经伸缩、时移、频移和频率倾斜的Gauss函数一一“chirplet”作为原子，以取代频率不变的Gabor原子(S.Mann另增加了时间倾斜操作)，并用内积法得到所谓的“chirplet变换”。

基于自适应匹配投影塔形分解算法的Chirplet变换，在本质上是对时频平面上的任意一条能量曲线用一组任意倾斜的线段进行线性逼近。显然，一阶逼近比零阶逼近能更紧凑地表达chirp类信号。尽管学者们一直热爱线性，然而大自然通常是非线性的。当信号的频率成分随时间非线性地变化时，譬如一类天然存在的或人工产生的Doppler信号，用频率随时间作线性(零阶或一阶)变化的原子来表征，势必造成原子数目的增加，从而既影响对分解结果的理解和途释，又影响分解结果的数据压缩能力。

但上述自适应时频分解方法基于已知的信号层次结构或者信号模型，在信号参数上是自适应分解的，但在信号类型上是非盲的。为进一步改进上述参数化时频分析方法还需要进行新的尝试。

时频分析在语音识别，雷达信号处理和图像处理，地震信号处理，信号重构以及扩频通信中的干扰抑制等方面，已有不少成功应用。总体说来，时频分析的应用领域大致有如下四类：一是时变谱分析；二是由时频分布间接计算出某些物理量；三是利用时频分布作为信号所携带的信息的载体(而不关心它是否真能表示能量密度)；四是信号的重构、压缩和编码等。

从考察信号的频率成分随时间的演化特性角度来说，小波变换的结果令人费解，尽管这一领域炙手可热。小波变换是以时间和尺度为参数，在时间-尺度平面的不同位置上具有不同的分辨率，因而是一种多分辨率分析方法。小波分析得益于小波原子的完备性、自相似性和多分辨性，它能获得成功的两个最重要的原因，是它拥有塔形快速算法和良好的时频局域特性；缺点则是一旦母小波选择不当，应用效果会大受影响。从信号压缩及消除交叉项干扰角度看，参数化的时频分析方法较好，但求其原子模型的参数，也非易事。WVD及所有其它的Cohen类时频分布，都可用于分析窄带信号，不过它们均不太适于分析多分量的宽带信号以及雷达和声纳信号，若采用仿射类的分布，则存在分辨率空变的问题。客观地说，各种时频分析技术难分优劣，关键是其适合何种类型的信号。一条屡试不爽的经验是，可以首先用运算速度很快的STFT的谱图尝试，如图3所示；若需要较高的时频分辨率，可以采用参数化时频分析方法。

但是，现有各种时频分析方法尚难完全适用于相位调制可表示为有限项多项式级数的多项式相位信号。

专利“一种模型驱动的多项式相位信号自适应时频变换方法”针对上述现实问题，提出了一种可行的新思路和新方法，但并未具体给出其中各种现代优化算法的具体的实用操作和运算过程。本发明针对这一问题，给出了一种用于非平稳多项式相位信号的自适应时频分解新算法，该算法是“模型驱动的多项式相位信号自适应时频变换方法”的一种具体实现方式和实用运算过程及计算流程，可以自适应的完成多项式相位信号的时频分解，并输出无任何交叉项干扰且时频联合分辨率较优的时频分布。

发明内容

为解决上述问题，本发明提供一种基于蚁群优化的多项式相位信号自适应时频变换方法，利用各信号分量、各时刻的瞬时频率取值，通过仅保留主瓣响应的Sinc函数直接计算生成相应时刻对应的信号频率分布，克服了传统的时频变换中一个时刻对应多个频点的非单分量存在交叉项的缺陷，最终输出无任何交叉项干扰且时频联合分辨率较优的时频分布。

一种基于蚁群优化的多项式相位信号自适应时频变换方法，包括以下步骤：

步骤1：利用蚁群优化算法对雷达接收的原始多项式相位信号s(t)进行多次分解，每次分解均确定一个信号分量的最优模型阶数和最优待定系数集，具体的：

步骤11：根据多项式相位模型，生成原始多项式相位信号s(t)对应的参考函数h_p，并初始化参考函数h_p的模型阶数N₁＝1，并计算原始多项式相位信号s(t)的能量E₀，随机生成参考函数h_p的至少两个待定系数集{a_n},n＝0,1,2,...,N₁，其中待定系数集{a_n}的个数用Ants表示；将待定系数集{a_n}作为蚁群优化算法中蚂蚁的位置参数PosAnts，Ants个蚂蚁构成初始蚂蚁群I_Ants；

步骤12：将各个待定系数集{a_n}对应的参考函数h_p的共轭分别与原始多项式相位信号s(t)相乘得到混合调制信号x(t)，对混合调制信号x(t)实施傅利叶变换，得到变换结果X(f)，计算变换结果X(f)其中一个频域特征的特征值，并将该特征值作为评价变换结果X(f)的适应度，根据各个变换结果X(f)的适应度优劣，得到当前优化周期的本地各个蚂蚁的位置参数AntL，本地各个蚂蚁的位置参数AntL对应的适应度Lfit，当前最优变换结果X(f)对应的全局最优蚂蚁的位置参数AntG，全局最优蚂蚁的位置参数AntG对应的适应度Gfit，其中AntG为所有优化周期中的最优蚂蚁的位置参数，且第一次优化周期中，AntL为初始蚂蚁群的位置参数PosAnts；

步骤13：根据设定的状态转移概率更新本地各个蚂蚁的位置参数，得到新蚂蚁群；

步骤14：采用所述新蚂蚁群中各蚂蚁位置参数对应的待定系数集{a_n}按步骤12的方法，重新计算并评价所述新蚂蚁群的位置参数对应的变换结果X(f)的适应度，并得到当前优化周期的本地各个蚂蚁的位置参数AntL、适应度Lfit、当前的全局最优蚂蚁的位置参数AntG以及适应度Gfit；

步骤15：对比步骤12中本地各个蚂蚁的位置参数的适应度Lfit和步骤14中新蚂蚁群各个蚂蚁的位置参数的适应度，如果更新后的适应度优于更新前的适应度，则本地各个蚂蚁采用更新后的位置参数，否则保留更新前的位置参数，对比完成后得到本地各个蚂蚁对应的适应度Part_fit；

步骤16：更新本地各个蚂蚁对应的适应度Part_fit；

步骤17：将更新后的本地各个蚂蚁对应的适应度Part_fit_new逐个代替步骤12中的适应度Gfit，重新得到状态转移概率，重复步骤13-步骤16，直至蚁群优化的周期达到设定的最大值Cycles，得到最终的全局最优蚂蚁的位置参数AntG；

步骤18：令模型阶数N₁依次从2取到最大可能阶数order_max，每次取值后，重复步骤11至步骤17，从而得到不同模型阶数下的经过Cycles个周期蚁群优化后的适应度最优的位置参数AntG；

步骤19：从步骤17得到所有位置参数AntG中选出适应度最优的位置参数，从而确定原始多项式相位信号s(t)当前信号分量的最优模型阶数N_p以及N_p对应的最优待定系数集{a_n}_max；

步骤2：利用最优模型阶数N_p和最优待定系数集{a_n}_max确定的信号分量的共轭与原始多项式相位信号s(t)相乘并实施傅利叶变换，得到频谱X'(f)，将频谱X'(f)包络最大值处的强度复数取值置零后，实施逆傅立叶变换，从而得到时域信号y(t)；

步骤3：利用模型阶数为N_p和最优待定系数集{a_n}_max确定的信号分量与时域信号y(t)相乘得到残差信号z(t)，本次分解结束；

步骤4：计算残差信号z(t)的能量E_d，与原始多项式相位信号s(t)能量E₀取比值R，如果比值R小于设定门限γ或分解次数达到设定的上限数量N_max，则停止分解，得到各分量信号的最优模型阶数N_p和最优待定系数集{a_n}_max；否则，利用步骤3计算的残差信号z(t)替换步骤1中的原始多项式相位信号s(t)重新计算混合调制信号x(t)，重复步骤1至步骤3，直至比值R小于设定门限γ或分解次数达到设定的上限数量N_max；

步骤5：将最优模型阶数N_p和最优待定系数集{a_n}_max对应的各个信号分量按分解顺序编号为h_c，其中c＝1,2,...,C，其中C为自适应分解的次数；

步骤6：按单分量信号瞬时频率的物理定义，得到各信号分量h_c对应的瞬时频率曲线f_c(t)，并根据所有信号分量h_c的瞬时频率曲线f_c(t)确定原始多项式相位信号s(t)的频率分布范围，最后对所述频率分布范围依次进行离散化、保留主瓣响应以及叠加操作后，得到原始多项式相位信号s(t)的时频联合分布f(t)。

进一步地，所述原始多项式相位信号s(t)对应的参考函数h_p具体为：

其中，order_max为原始多项式相位信号s(t)的最大可能阶数。

进一步地，所述原始多项式相位信号s(t)的频率分布范围具体为：

所有信号分量h_c的瞬时频率曲线f_c(t)中的瞬时频率最小值为原始多项式相位信号s(t)的最小频率f_imin，瞬时频率最大值为原始多项式相位信号s(t)的最大频率f_imax。

进一步地，所述待定系数集{a_n}具体为矢量

其中N₁代表当前优化的模型阶数，a₀为固定取值，设为原始多项式相位信号s(t)直接傅利叶变换所得峰值的所在采样点的位置序号；

为待优化的待定系数，由随机函数Rand产生，取值范围为[-M,M]，其中M为原始多项式相位信号s(t)采样序列的采样点数目，具体的：

其中rand(Ants,N₁)表示随机产生Ants个维数为N₁的分布于[0,1]的随机数，从而得到模型阶数为N₁的初始蚂蚁群I_Ants。

进一步地，步骤12所述的频域特征为频谱峰值、频谱的对比度或频谱的谱形熵中的任意一个，其中各频域特征的特征值频谱峰值X_p、频谱的对比度X_c以及频谱的谱形熵X_e分别为：

X_p＝max{abs[X(f)]}

X_e＝sum(-p_S·logp_S)

其中

sum(·)代表向量的求和操作，abs(·)为取绝对值；

当选择频谱峰值X_p或频谱对比度X_c作为适应度时，该值越大，评价结果越优，当选择频谱的谱形熵X_e作为适应度时，该值越小，评价结果越优。

进一步地，步骤13所述的根据设定的状态转移概率更新步骤12中本地各个蚂蚁的位置参数，具体计算方法为：

P_ant＝(Gfit-Lfit)/Gfit

其中，P_ant为状态转移概率，P_ant0为设定的转移门限，P_step为蚁群优化设定的局部搜索步长，λ为蚁群优化周期次数的倒数，具体的：

其中cyc取值范围为1,2,...,Cycles。

进一步地，对步骤13更新后新蚂蚁群的位置参数进行纠错，完成纠错的位置参数再进入步骤14，其中纠错具体采用的规则为状态转移规则或边界规则；

当采用状态转移规则对更新后新蚂蚁群的位置参数进行纠错时，如果新蚂蚁群中任一蚂蚁其中一维的位置参数超过设定的取值范围[-M,M]，该蚂蚁的位置状态转移概率立即置为1，且该蚂蚁该维的位置参数将在下一周期优化过程中重新更新后返回求解范围；

当采用边界规则对更新后新蚂蚁群的位置参数进行纠错时，如果新蚂蚁群中任一蚂蚁其中一维的位置参数超过设定的取值范围[-M,M]，该蚂蚁该维的位置参数根据符号被强制置为位置边界取值M或-M，且赋予该蚂蚁随机的状态转移概率，同时该蚂蚁该维的位置参数将在下一周期优化过程中重新更新后返回求解范围；其中，随机的状态转移概率赋予规则为：P_ant＝rand，其中rand表示只产生一个分布在[0,1]范围内的随机数。

进一步地，步骤16所述的更新本地各个蚂蚁对应的适应度Part_fit，具体为根据设定的蚁群优化的信息素蒸发概率R_ho和步骤12中本地各个蚂蚁的位置参数AntL对应的适应度Lfit，实现对本地各个蚂蚁对应的适应度Part_fit的更新，具体的：

Part_fit_new＝(1-R_ho)×Lfit+Part_fit。

进一步地，步骤6所述的对所述频率分布范围依次进行离散化具体为：

f_is＝[f_imin f_imin+Δf f_imin+2Δf f_imin+3Δf ... f_imax]^T

其中，f_is为离散频率向量的瞬时频率值，Δf为设定的频率分辨率，具体的：

其中，f_imin为原始多项式相位信号s(t)的最小频率，f_imax为原始多项式相位信号s(t)的最大频率，K为设定的离散频率维数，i为虚部单位。

进一步地，步骤6所述的对所述频率分布范围进行保留主瓣响应具体计算方法为：

其中，A_c为第c个信号分量h_c的频谱包络最大值处的强度复数取值，i为虚部单位，Δf为设定的频率分辨率，f_is为离散频率向量f_s的瞬时频率值、f_ic(t)为各信号分量h_c对应的瞬时频率曲线f_c(t)的瞬时频率值。

有益效果：

本发明提供一种基于蚁群优化的多项式相位信号自适应时频变换方法，能够完成多项式相位信号的时频分解，其中分解得到的每一个信号分量都为任一时刻都只对应一个频点的单分量，然后利用各信号分量、各时刻的瞬时频率取值，通过仅保留主瓣响应的Sinc函数直接计算生成相应时刻对应的信号频率分布，克服了传统的时频变换中一个时刻对应多个频点的非单分量存在交叉项的缺陷，最终输出无任何交叉项干扰且时频联合分辨率较优的时频分布；

本发明原理简单，操作方便，可有效克服经典时频分析方法交叉项干扰的不利影响以及时频联合分辨率的损失，能够有效提升非平稳多项式相位信号时频分析的质量和效益。

附图说明

图1为本发明蚁群优化的自适应时频变换方法流程图；

图2为本发明其中一个多分量多项式相位信号的时域波形图；

图3为现有技术中经典STFT分析时频图；

图4为现有技术中经典WVD时频图；

图5为本发明蚁群优化的自适应时频变换方法所得时频图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细叙述。

本发明首先通过一定周期且具有一定规模的蚁群优化估计出多项式相位信号各组成分量的模型阶数及对应的各阶相位参数，然后利用“清洁”的思想，提取出相应的多项式相位信号分量，并从原始信号中剔除该信号分量得到残差信号，接着以迭代“清洁”的方式，重复利用蚁群优化针对残差信号实施自适应模型定阶及最优模型参数优化，并逐步提取出各信号分量，如此反复，直至残差信号能量低于预设的门限或提取出的信号分量数目达到预设的最大值。随后，利用蚁群优化分解所得的各分量信号的相位参数构建相位-时间历史，并根据单分量信号瞬时频率的物理定义，对相位历史求导，以解析的方式直接得到各信号分量的频率-时间历史，即时频变化曲线，根据信号各分量的时频变化曲线，确定整个信号的最大及最小频率，利用其确定整个信号的频率变化范围，按应用所需的频率分辨率需求，对该频率范围进行离散化，得到所需维数的离散频率向量，最后，利用各分量、各时刻的瞬时频率取值，通过仅保留主瓣响应的Sinc函数直接计算生成相应时刻对应的信号频率分布；如此反复，直至所有信号分量所有时刻的频率分布均生成完毕，将其按时间先后顺序存储并逐信号分量累加，即可得到最终多项式相位信号的时频联合分布。

如图1所示，一种基于蚁群优化的多项式相位信号自适应时频变换方法，包括以下步骤：

步骤1：利用蚁群优化算法对原始多项式相位信号s(t)进行多次分解，每次分解均确定一个信号分量的最优模型阶数和最优待定系数集，具体的：

步骤11：根据多项式相位模型，生成原始多项式相位信号s(t)对应的参考函数h_p，并初始化参考函数h_p的模型阶数N₁＝1，并计算原始多项式相位信号s(t)的能量E₀，随机生成参考函数h_p的至少两个待定系数集{a_n},n＝0,1,2,...,N₁，其中待定系数集{a_n}的个数用Ants表示；将待定系数集{a_n}作为蚁群优化算法中蚂蚁的位置参数PosAnts，Ants个蚂蚁构成初始蚂蚁群I_Ants；

步骤12：将各个待定系数集{a_n}对应的参考函数h_p的共轭分别与原始多项式相位信号s(t)相乘得到混合调制信号x(t)，对混合调制信号x(t)实施傅利叶变换，得到变换结果X(f)，计算变换结果X(f)其中一个频域特征的特征值，并将该特征值作为评价变换结果X(f)的适应度，根据各个变换结果X(f)的适应度优劣，得到当前优化周期的本地各个蚂蚁的位置参数AntL，本地各个蚂蚁的位置参数AntL对应的适应度Lfit，当前最优变换结果X(f)对应的全局最优蚂蚁的位置参数AntG，全局最优蚂蚁的位置参数AntG对应的适应度Gfit，其中AntG为所有优化周期中的最优蚂蚁的位置参数，且第一次优化周期中，AntL为初始蚂蚁群的位置参数PosAnts；

步骤13：根据设定的状态转移概率更新本地各个蚂蚁的位置参数，得到新蚂蚁群；

步骤14：采用所述新蚂蚁群中各蚂蚁位置参数对应的待定系数集{a_n}按步骤12的方法，重新计算并评价所述新蚂蚁群的位置参数对应的变换结果X(f)的适应度，并得到当前优化周期的本地各个蚂蚁的位置参数AntL、适应度Lfit、当前的全局最优蚂蚁的位置参数AntG以及适应度Gfit；

步骤15：对比步骤12中本地各个蚂蚁的位置参数的适应度Lfit和步骤14中新蚂蚁群各个蚂蚁的位置参数的适应度，如果更新后的适应度优于更新前的适应度，则本地各个蚂蚁采用更新后的位置参数，否则保留更新前的位置参数，对比完成后得到本地各个蚂蚁对应的适应度Part_fit；

步骤16：根据设定的蚁群优化的信息素蒸发概率R_ho和步骤12中本地各个蚂蚁的位置参数AntL对应的适应度Lfit，更新本地各个蚂蚁对应的适应度Part_fit；

步骤17：将更新后的本地各个蚂蚁对应的适应度Part_fit_new分别代替步骤12中的适应度Gfit，重新得到状态转移概率，重复步骤13-步骤16，直至蚁群优化的周期达到设定的最大值Cycles，得到最终的全局最优蚂蚁的位置参数AntG；

步骤18：令模型阶数N₁依次从2取到最大可能阶数order_max，每次取值后，重复步骤11至步骤17，从而得到不同模型阶数下的经过Cycles个周期蚁群优化后的适应度最优的位置参数AntG；

步骤19：从步骤17得到所有位置参数AntG中选出适应度最优的位置参数，从而确定原始多项式相位信号s(t)当前信号分量的最优模型阶数N_p以及N_p对应的最优待定系数集{a_n}_max；

步骤2：利用最优模型阶数N_p和最优待定系数集{a_n}_max确定的信号分量的共轭与原始多项式相位信号s(t)相乘并实施傅利叶变换，得到频谱X'(f)，将频谱X'(f)包络最大值处的强度复数取值置零后，实施逆傅立叶变换，从而得到时域信号y(t)，如图2所示；

步骤3：利用模型阶数为N_p和最优待定系数集{a_n}_max确定的信号分量与时域信号y(t)相乘得到残差信号z(t)，本次分解结束；

步骤4：计算残差信号z(t)的能量E_d，与原始多项式相位信号s(t)能量E₀取比值R，如果比值R小于设定门限γ或分解次数达到设定的上限数量N_max，则停止分解，得到各分量信号的最优模型阶数N_p和最优待定系数集{a_n}_max；否则，利用步骤3计算的残差信号z(t)替换步骤1中的原始多项式相位信号s(t)重新计算混合调制信号x(t)，重复步骤1至步骤3，直至比值R小于设定门限γ或分解次数达到设定的上限数量N_max；

步骤5：将最优模型阶数N_p和最优待定系数集{a_n}_max对应的各个信号分量按分解顺序编号为h_c，其中c＝1,2,...,C，其中C为自适应分解的次数；

步骤6：按单分量信号瞬时频率的物理定义，得到各信号分量h_c对应的瞬时频率曲线f_c(t)，并根据所有信号分量h_c的瞬时频率曲线f_c(t)确定原始多项式相位信号s(t)的频率分布范围，最后对所述频率分布范围依次进行离散化、保留主瓣响应以及叠加操作后，得到原始多项式相位信号s(t)的时频联合分布f(t)，如图5所示。

进一步地，所述原始多项式相位信号s(t)对应的参考函数h_p具体为：

其中，待定系数集{a_n},n＝0,1,2,...N₁,中的待定系数a_n初始值随机生成，order_max为原始多项式相位信号s(t)的最大可能阶数。

进一步地，所述原始多项式相位信号s(t)的频率分布范围具体为：

所有信号分量h_c的瞬时频率曲线f_c(t)中的瞬时频率最小值为原始多项式相位信号s(t)的最小频率f_{i min}，瞬时频率最大值为原始多项式相位信号s(t)的最大频率f_{i max}。

进一步地，所述待定系数集{a_n}具体为矢量

其中N₁代表当前优化的模型阶数，a₀为固定取值，设为原始多项式相位信号s(t)直接傅利叶变换所得峰值的所在采样点的位置序号；

为待优化的待定系数，由随机函数Rand产生，取值范围为[-M,M]，其中M为原始多项式相位信号s(t)采样序列的采样点数目，具体的：

其中rand(Ants,N₁)表示随机产生Ants个维数为N₁的分布于[0,1]的随机数，从而得到模型阶数为N₁的初始蚂蚁群I_Ants。

进一步地，步骤12所述的频域特征为频谱峰值、频谱的对比度或频谱的谱形熵中的任意一个，其中各频域特征的特征值频谱峰值X_p、频谱的对比度X_c以及频谱的谱形熵X_e分别为：

X_p＝max{abs[X(f)]}

X_e＝sum(-p_S·logp_S)

其中

sum(·)代表向量的求和操作，abs(·)为取绝对值；

当选择频谱峰值X_p或频谱对比度X_c作为适应度时，该值越大，评价结果越优，当选择频谱的谱形熵X_e作为适应度时，该值越小，评价结果越优。

进一步地，步骤13所述根据设定的状态转移概率更新步骤12中本地各个蚂蚁的位置参数，具体计算方法为：

P_ant＝(Gfit-Lfit)/Gfit

其中，P_ant为状态转移概率，P_ant0为设定的转移门限，P_step为蚁群优化设定的局部搜索步长，λ为蚁群优化周期次数的倒数，具体的：

其中cyc取值范围为1,2,...,Cycles。

进一步地，对步骤13更新后新蚂蚁群的位置参数进行纠错，完成纠错的位置参数再进入步骤14，其中纠错具体采用的规则为状态转移规则或边界规则；

当采用状态转移规则对更新后新蚂蚁群的位置参数进行纠错时，如果新蚂蚁群中任一蚂蚁其中一维的位置参数超过设定的取值范围[-M,M]，该蚂蚁的位置状态转移概率立即置为1，且该蚂蚁该维的位置参数将在下一周期优化过程中重新更新后返回求解范围；

当采用边界规则对更新后新蚂蚁群的位置参数进行纠错时，如果新蚂蚁群中任一蚂蚁其中一维的位置参数超过设定的取值范围[-M,M]，该蚂蚁该维的位置参数根据符号被强制置为位置边界取值M或-M，且赋予该蚂蚁随机的状态转移概率，同时该蚂蚁该维的位置参数将在下一周期优化过程中重新更新后返回求解范围；其中，随机的状态转移概率赋予规则为：P_ant＝rand，其中rand表示只产生一个分布在[0,1]范围内的随机数。

进一步地，步骤16所述的更新本地各个蚂蚁对应的适应度Part_fit，具体为：

Part_fit_new＝(1-R_ho)·Lfit+Part_fit。

进一步地，步骤6所述的对所述频率分布范围依次进行离散化具体为：

f_is＝[f_imin f_imin+Δf f_imin+2Δf f_imin+3Δf ... f_imax]^T

其中，f_is为离散频率向量的瞬时频率值，Δf为设定的频率分辨率，具体的：

其中，f_imin为原始多项式相位信号s(t)的最小频率，f_imax为原始多项式相位信号s(t)的最大频率，K为设定的离散频率维数，i为虚部单位。

进一步地，步骤6所述的对所述频率分布范围进行保留主瓣响应具体计算方法为：

其中，A_c为第c个信号分量h_c的频谱包络最大值处的强度复数取值，i为虚部单位，Δf为设定的频率分辨率，f_is为离散频率向量f_s的瞬时频率值、f_ic(t)为各信号分量h_c对应的瞬时频率曲线f_c(t)的瞬时频率值。

当然，本发明还可有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (10)

1.一种基于蚁群优化的多项式相位信号自适应时频变换方法，其特征在于，包括以下步骤：

步骤1：利用蚁群优化算法对雷达接收的原始多项式相位信号s(t)进行多次分解，每次分解均确定一个信号分量的最优模型阶数和最优待定系数集，具体的：

步骤11：根据多项式相位模型，生成原始多项式相位信号s(t)对应的参考函数h_p，并初始化参考函数h_p的模型阶数N₁＝1，并计算原始多项式相位信号s(t)的能量E₀，随机生成参考函数h_p的至少两个待定系数集{a_n},n＝0,1,2,...,N₁，其中待定系数集{a_n}的个数用Ants表示；将待定系数集{a_n}作为蚁群优化算法中蚂蚁的位置参数PosAnts，Ants个蚂蚁构成初始蚂蚁群I_Ants；

步骤12：将各个待定系数集{a_n}对应的参考函数h_p的共轭分别与原始多项式相位信号s(t)相乘得到混合调制信号x(t)，对混合调制信号x(t)实施傅利叶变换，得到变换结果X(f)，计算变换结果X(f)其中一个频域特征的特征值，并将该特征值作为评价变换结果X(f)的适应度，根据各个变换结果X(f)的适应度优劣，得到当前优化周期的本地各个蚂蚁的位置参数AntL，本地各个蚂蚁的位置参数AntL对应的适应度Lfit，当前最优变换结果X(f)对应的全局最优蚂蚁的位置参数AntG，全局最优蚂蚁的位置参数AntG对应的适应度Gfit，其中AntG为所有优化周期中的最优蚂蚁的位置参数，且第一次优化周期中，AntL为初始蚂蚁群的位置参数PosAnts；

步骤13：根据设定的状态转移概率更新本地各个蚂蚁的位置参数，得到新蚂蚁群；

步骤14：采用所述新蚂蚁群中各蚂蚁位置参数对应的待定系数集{a_n}按步骤12的方法，重新计算并评价所述新蚂蚁群的位置参数对应的变换结果X(f)的适应度，并得到当前优化周期的本地各个蚂蚁的位置参数AntL、适应度Lfit、当前的全局最优蚂蚁的位置参数AntG以及适应度Gfit；

步骤15：对比步骤12中本地各个蚂蚁的位置参数的适应度Lfit和步骤14中新蚂蚁群各个蚂蚁的位置参数的适应度，如果更新后的适应度优于更新前的适应度，则本地各个蚂蚁采用更新后的位置参数，否则保留更新前的位置参数，对比完成后得到本地各个蚂蚁对应的适应度Part_fit；

步骤16：更新本地各个蚂蚁对应的适应度Part_fit；

步骤17：将更新后的本地各个蚂蚁对应的适应度Part_fit_new逐个代替步骤12中的适应度Gfit，重新得到状态转移概率，重复步骤13-步骤16，直至蚁群优化的周期达到设定的最大值Cycles，得到最终的全局最优蚂蚁的位置参数AntG；

步骤18：令模型阶数N₁依次从2取到最大可能阶数order_max，每次取值后，重复步骤11至步骤17，从而得到不同模型阶数下的经过Cycles个周期蚁群优化后的适应度最优的位置参数AntG；

步骤19：从步骤17得到所有位置参数AntG中选出适应度最优的位置参数，从而确定原始多项式相位信号s(t)当前信号分量的最优模型阶数N_p以及N_p对应的最优待定系数集{a_n}_max；

步骤2：利用最优模型阶数N_p和最优待定系数集{a_n}_max确定的信号分量的共轭与原始多项式相位信号s(t)相乘并实施傅利叶变换，得到频谱X'(f)，将频谱X'(f)包络最大值处的强度复数取值置零后，实施逆傅利叶变换，从而得到时域信号y(t)；

步骤3：利用模型阶数为N_p和最优待定系数集{a_n}_max确定的信号分量与时域信号y(t)相乘得到残差信号z(t)，本次分解结束；

步骤4：计算残差信号z(t)的能量E_d，与原始多项式相位信号s(t)能量E₀取比值R，如果比值R小于设定门限γ或分解次数达到设定的上限数量N_max，则停止分解，得到各分量信号的最优模型阶数N_p和最优待定系数集{a_n}_max；否则，利用步骤3计算的残差信号z(t)替换步骤1中的原始多项式相位信号s(t)重新计算混合调制信号x(t)，重复步骤1至步骤3，直至比值R小于设定门限γ或分解次数达到设定的上限数量N_max；

步骤5：将最优模型阶数N_p和最优待定系数集{a_n}_max对应的各个信号分量按分解顺序编号为h_c，其中c＝1,2,...,C，其中C为自适应分解的次数；

步骤6：按单分量信号瞬时频率的物理定义，得到各信号分量h_c对应的瞬时频率曲线f_c(t)，并根据所有信号分量h_c的瞬时频率曲线f_c(t)确定原始多项式相位信号s(t)的频率分布范围，最后对所述频率分布范围依次进行离散化、保留主瓣响应以及叠加操作后，得到原始多项式相位信号s(t)的时频联合分布f(t)。

2.如权利要求1所述的一种基于蚁群优化的多项式相位信号自适应时频变换方法，其特征在于，所述原始多项式相位信号s(t)对应的参考函数h_p具体为：

其中，order_max为原始多项式相位信号s(t)的最大可能阶数，j为虚部单位。

3.如权利要求1所述的一种基于蚁群优化的多项式相位信号自适应时频变换方法，其特征在于，所述原始多项式相位信号s(t)的频率分布范围具体为：

所有信号分量h_c的瞬时频率曲线f_c(t)中的瞬时频率最小值为原始多项式相位信号s(t)的最小频率f_imin，瞬时频率最大值为原始多项式相位信号s(t)的最大频率f_imax。

4.如权利要求1所述的一种基于蚁群优化的多项式相位信号自适应时频变换方法，其特征在于，所述待定系数集{a_n}具体为矢量

其中N₁代表当前优化的模型阶数，a₀为固定取值，设为原始多项式相位信号s(t)直接傅利叶变换所得峰值的所在采样点的位置序号；

为待优化的待定系数，由随机函数Rand产生，取值范围为[-M,M]，其中M为原始多项式相位信号s(t)采样序列的采样点数目，具体的：

其中rand(Ants,N₁)表示随机产生Ants个维数为N₁的分布于[0,1]的随机数，从而得到模型阶数为N₁的初始蚂蚁群I_Ants。

5.如权利要求1所述的一种基于蚁群优化的多项式相位信号自适应时频变换方法，其特征在于，步骤12所述的频域特征为频谱峰值、频谱的对比度或频谱的谱形熵中的任意一个，其中各频域特征的特征值频谱峰值X_p、频谱的对比度X_c以及频谱的谱形熵X_e分别为：

X_p＝max{abs[X(f)]}

X_e＝sum(-p_S·logp_S)

其中

sum(·)代表向量的求和操作，abs(·)为取绝对值；

当选择频谱峰值X_p或频谱对比度X_c作为适应度时，该值越大，评价结果越优，当选择频谱的谱形熵X_e作为适应度时，该值越小，评价结果越优。

6.如权利要求1所述的一种基于蚁群优化的多项式相位信号自适应时频变换方法，其特征在于，步骤13所述的根据设定的状态转移概率更新步骤12中本地各个蚂蚁的位置参数，具体计算方法为：

P_ant＝(Gfit-Lfit)/Gfit

其中，P_ant为状态转移概率，P_ant0为设定的转移门限，P_step为蚁群优化设定的局部搜索步长，λ为蚁群优化周期次数的倒数，具体的：

其中cyc取值范围为1,2,...,Cycles。

7.如权利要求1所述的一种基于蚁群优化的多项式相位信号自适应时频变换方法，其特征在于，对步骤13更新后新蚂蚁群的位置参数进行纠错，完成纠错的位置参数再进入步骤14，其中纠错具体采用的规则为状态转移规则或边界规则；

当采用状态转移规则对更新后新蚂蚁群的位置参数进行纠错时，如果新蚂蚁群中任一蚂蚁其中一维的位置参数超过设定的取值范围[-M,M]，该蚂蚁的位置状态转移概率立即置为1，且该蚂蚁该维的位置参数将在下一周期优化过程中重新更新后返回求解范围；

当采用边界规则对更新后新蚂蚁群的位置参数进行纠错时，如果新蚂蚁群中任一蚂蚁其中一维的位置参数超过设定的取值范围[-M,M]，该蚂蚁该维的位置参数根据符号被强制置为位置边界取值M或-M，且赋予该蚂蚁随机的状态转移概率，同时该蚂蚁该维的位置参数将在下一周期优化过程中重新更新后返回求解范围；其中，随机的状态转移概率赋予规则为：P_ant＝rand，其中rand表示只产生一个分布在[0,1]范围内的随机数。

8.如权利要求1所述的一种基于蚁群优化的多项式相位信号自适应时频变换方法，其特征在于，步骤16所述的更新本地各个蚂蚁对应的适应度Part_fit，具体为根据设定的蚁群优化的信息素蒸发概率R_ho和步骤12中本地各个蚂蚁的位置参数AntL对应的适应度Lfit，实现对本地各个蚂蚁对应的适应度Part_fit的更新，具体的：

Part_fit_new＝(1-R_ho)×Lfit+Part_fit。

9.如权利要求1所述的一种基于蚁群优化的多项式相位信号自适应时频变换方法，其特征在于，步骤6所述的对所述频率分布范围依次进行离散化具体为：

f_is＝[f_imin f_imin+Δf f_imin+2Δf f_imin+3Δf...f_imax]^T

其中，f_is为离散频率向量的瞬时频率值，Δf为设定的频率分辨率，具体的：

其中，f_imin为原始多项式相位信号s(t)的最小频率，f_imax为原始多项式相位信号s(t)的最大频率，K为设定的离散频率维数，i为虚部单位。

10.如权利要求1所述的一种基于蚁群优化的多项式相位信号自适应时频变换方法，其特征在于，步骤6所述的对所述频率分布范围进行保留主瓣响应具体计算方法为：

其中，A_c为第c个信号分量h_c的频谱包络最大值处的强度复数取值，i为虚部单位，Δf为设定的频率分辨率，f_is为离散频率向量f_s的瞬时频率值、f_ic(t)为各信号分量h_c对应的瞬时频率曲线f_c(t)的瞬时频率值。

Images