1.一种基于前-后向代价参考粒子滤波的瞬时频率曲线估计方法，其特征在于，包括以下步骤：
(1)将观测时间区间[0,T]等分为L个不相交的子区间：[(k-1)ΔT,kΔT]，构建第k个观测时刻信号状态向量xk的状态方程和观测向量zk的观测方程：
其中，ΔT＝T/L，A是状态转移矩阵，E噪声转移矩阵，vk表示第k时刻统计特性未知的噪声向量，ak表示第k时刻回波信号的幅度， 表示第k时刻回波信号的初始相位，h(xk)表示幅度为1的线性调时频粒子，wk表示第k时刻的观测噪声，该测噪声服从均值为零，方
2
差为σ的高斯分布；
(2)用户根据状态向量xk和观测向量zk定义粒子代价ck：
2
其中 表示二范数的平方，|·|表示模平方， 表示zk的共轭转置，M表示观测向量的长度，ck的值越小，则表明h(xk)近似观测zk的程度越好；
(3)将第0个观测时刻作为初始时刻，在状态向量xk所属的状态空间Ω中均匀采样产生N个初始粒子 设每个初始粒子的初始代价为0，即 得到
初始时刻的粒子-代价集合为
(4)根据第k时刻的粒子-代价集合 通过步骤(1)中的状
态方程和观测方程预测出第k+1时刻的粒子 计算k时刻第n个粒子的风险
其中，0≤λ＜1为遗忘因子，zk+1为k+1时刻的观测向量， 值越小的粒子对下一时刻状态的估计越重要；
(5)根据步骤(4)的粒子风险 计算重采样权值：
其中q≥1是一个经验参数；
(6)在k+1时刻，根据k时刻N个粒子的重采样权值之和进行重采样和粒子代价更新，得到k+1时刻的粒子-代价集合Λk+1；
(7)重复步骤(4)至步骤(6)直至迭代到最后时刻L，得到L个观测时刻的粒子-代价集合：Λ1,Λ2,…,ΛL，以每个时刻最小的代价为估计准则，得到基于前向代价参考粒子滤波的信号状态估计
其中 表示 最小时n的取值，n0表示代价最小时对应粒子的序号；
(8)构建后向动态系统信号的状态方程和观测方程：
其中，AB为后向动态系统的状态转移矩阵，EB为噪声转移矩阵，vk-1表示第k-1时刻统计特性未知的噪声向量；
(9)根据步骤(8)确定的后向动态系统信号的状态方程和观测方程，将步骤(7)得到的第L时刻的粒子-代价集合作为初始粒子-代价集合，并将观测数据逆序代入后向动态系统的状态方程和观测方程，得到L个前-后向动态系统的粒子-代价集合以每个时刻最小的代价为估计准则，得到基于前-后向代价参考粒子
滤波的信号状态估计
其中， 表示前-后向动态系统中k时刻的
第n个粒子， 表示粒子 的代价， 表示中心频率的估计， 表示频率变化率的估计；
(10)利用基于前-后向代价参考粒子滤波的信号状态估计 中的中心频率估计 与频率变化率估计 得到信号瞬时频率曲线的估计
2.如权利要求1所述的方法，其中所述步骤(1)中构建第k个观测时刻信号状态向量xk的状态方程和观测向量zk的观测方程，按如下步骤进行：
1a)将观测时间区间[0,T]等分为L个不相交的子区间：[(k-1)ΔT,kΔT]，在第k个观测时刻，用线性调频时频粒子h(xk)模拟非线性调频信号sk(t)：
其中k＝1,2,...,L，ak表示第k时刻信号的幅度， 表示第k时刻信号的初始相位，fk表示第k时刻信号的中心频率，rk表示第k时刻信号的频率变化率；
1b)用幅度ak、初始相位 中心频率fk和频率变化率rk四个参数来描述每个子区间中信号sk(t)的状态，并忽略幅度ak和初始相位 得到第k时刻信号的状态向量：xk＝[fk,rk]；
1c)将信号瞬时频率曲线fd(t)及其导数f′d(t)在tk-1＝(k-1.5)ΔT处进行泰勒展开，得到展开后k时刻的中心频率fk和频率变化率rk：
rk＝f′ d((k-0.5)ΔT)＝rk-1+ΔTfd″(ξ2)
其中，fd((k-0.5)ΔT)表示fd(t)在tk＝(k-0.5)ΔT处的值，f′ d((k-0.5)ΔT)表示f′d(t)在tk＝(k-0.5)ΔT处的值，ξ1、ξ2∈((k-1.5)ΔT,(k-0.5)ΔT)，ΔTfd″(ξ1)是fd(t)在tk-1处泰勒展开式的余项，ΔTfd″(ξ2)是f′d(t)在tk-1处泰勒展开式的余项，ΔTfd″(ξ1)与ΔTfd″(ξ2)表示两个统计特性未知的噪声变量；
1d)根据步骤1c)泰勒展开后的中心频率fk和频率变化率rk，得到信号状态向量xk的状态方程和观测向量zk的观测方程：
其中， 表示第k时刻统计特性未知的噪
声向量。
3.如权利要求1所述的方法，其中所述步骤(6)中根据k时刻N个粒子的重采样权值之和进行重采样和粒子代价更新，按如下步骤进行：
6a)计算k时刻N个粒子的重采样权值之和： 其中，n＝1,2,…,N，[·]
表示取整运算；
6b)在k+1时刻，根据αn进行重采样和粒子代价更新：
当αn≥1时，利用粒子 重采样产生αn个重采样粒子： l＝
1,2,…,αn，并更新每个重采样粒子的代价为：
其中， 表示重采样粒子， 表示重采样粒子的代价， 是随机数；
当αn＜ 1时，不进行重采样，则 k+1时刻的粒子 -代价集合为：
4.如权利要求1所述的方法，其中所述步骤(8)中构建后向动态系统信号的状态方程和观测方程，按如下步骤进行：
8a)将信号瞬时频率曲线fd(t)及其导数f′d(t)在tk＝(k-0.5)ΔT处进行泰勒展开，得到展开后k-1时刻的中心频率fk-1和频率变化率rk-1：
2
fk-1＝fd((k-1.5)ΔT)＝fk-ΔTrk-1+0.5(ΔT)fd″(ξ3)，
rk-1＝f′ d((k-1.5)ΔT)＝rk-ΔTfd″(ξ4)
其中，fd((k-1.5)ΔT)表示fd(t)在tk＝(k-1.5)ΔT处的值，f′ d((k-1.5)ΔT)表示f′d(t)在tk＝(k-1.5)ΔT处的值，ξ3、ξ4∈((k-1.5)ΔT,(k-0.5)ΔT)，ΔTfd″(ξ3)是fd(t)在tk处泰勒展开式的余项，ΔTfd″(ξ4)是f′d(t)在tk处泰勒展开式的余项，ΔTfd″(ξ3)与ΔTfd″(ξ4)表示两个统计特性未知的噪声变量；
8b)根据步骤8a)泰勒展开后的中心频率fk-1和频率变化率rk-1，得到后向动态系统信号的状态方程和观测方程：
其中 表示k-1时刻统计特性未知的
噪声向量。