CN101894097B - 卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法 - Google Patents

Abstract

卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法，涉及机动目标滤波跟踪领域，解决了无法克服EMD方法中的“边界效应”及强机动目标卡尔曼滤波的发散问题。该方法为：一、获得当前采样时刻的量测数据序列，利用卡尔曼预测方程，获得下一时刻的系统状态预测数据，将其与当前序列组合；二、在组合序列中插值生成滤波序列；三、利用EMD方法对滤波序列进行分解，剔除包含噪声的IMF，获得当前滤波值；四、将该滤波值作为当前滤波结果显示；五、将该滤波值作为系统当前状态的后验估计，在下一时刻结合量测数据及卡尔曼方程组，获得计算参数，然后返回一，将下一时刻作为当前时刻，实现机动目标的跟踪。本发明可用于机动目标滤波跟踪领域。

Description

卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法

技术领域

本发明涉及机动目标滤波跟踪问题，具体涉及一种卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法。

背景技术

目标跟踪滤波方法在探测、识别与跟踪等多个领域有着重要现实意义和广阔的应用背景，其中机动目标跟踪滤波问题更是研究热点之一。经过众多学者的不懈努力提出或发展了许多相关的理论，如卡尔曼理论及其扩展方法、小波分析、粒子滤波、多模型概率估计等等，其中以卡尔曼理论及其扩展的滤波方法最为人们接收，应用也是最为广泛，但它们和许多其他方法一样，仍存在着对参数依赖性强、易发散等问题，从一定程度上限制了这些方法的应用。

经验模态分解(Empirical Mode Decomposition-EMD)方法是于1998年由美国国家工程学院院士黄锷提出，该方法是将原信号分解为从高频到低频的固有模态函数(IntrinsicMode Function-IMF)，且分解出的各个IMF的频率高低是由信号本身的频率特性决定的，具有自适应的特性，因此适用于噪声形式未知的非平稳信号滤波。在应用EMD方法时，有一个无法回避的问题，即在用极点拟合上下包络的时候，在两个端点附近就会出现包络失真问题，造成计算上的不准确，这种现象被称为“边界效应”。而且随着求取IMF阶次的升高，迭代的次数不断增加，“边界效应”还会传播到信号内部，进一步影响EMD分解的质量，使分解出的IMF失去了实际的物理意义。而在对目标进行跟踪滤波时所关心的是目标最新的状态，恰恰是信号序列的边界值，正是这个原因限制了EMD在机动目标跟踪滤波中的应用。

发明内容

本发明的目的是克服EMD方法中的“边界效应”以及强机动目标卡尔曼滤波的发散问题，提供了一种卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法。

卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法，它的具体过程如下：

步骤一、获得当前采样时刻的量测数据序列，并利用卡尔曼预测方程，获得下一采样时刻的系统状态预测数据，并将所述下一采样时刻的系统状态预测数据和当前量测数据序列组合，获得组合后的数据序列；

步骤二、在步骤一获得的组合后的数据序列中的每相邻两个数据之间均匀分布多个插值点，生成滤波序列；

步骤三、利用EMD方法对步骤二生成的滤波序列进行分解，获得包含噪声的IMF，在滤波序列中剔除包含噪声的IMF，然后获得当前时刻的滤波值；

步骤四、将当前时刻的滤波值作为当前滤波结果显示；

步骤五、将当前时刻的滤波值作为系统当前状态的后验估计，在下一采样时刻来临时，结合下一采样时刻的量测数据，并利用卡尔曼方程组得到用于下一采样时刻的计算参数，然后，返回执行步骤一，将下一采样时刻作为当前采样时刻，实现机动目标的跟踪。

本发明的积极效果：

本发明有效地克服了EMD方法中“边界效应”以及针对强机动目标卡尔曼滤波的发散问题，具有较高的跟踪精度。

附图说明

图1为本发明的流程图；图2为步骤三的流程图；图3为步骤三一的流程图；图4为应用EKF方法在30秒时间段内测得的地理坐标系X轴向的目标速度信息的滤波效果图；图5为应用KF方法在30秒时间段内测得的地理坐标系X轴向的目标速度信息的滤波效果图；图6为应用EMD-KF方法在30秒时间段内测得的地理坐标系X轴向的目标速度信息的滤波效果图；图7为应用EKF方法在30秒时间段内测得的地理坐标系Y轴向的目标速度信息的滤波效果图；图8为应用KF方法在30秒时间段内测得的地理坐标系Y轴向的目标速度信息的滤波效果图；图9为应用EMD-KF方法在30秒时间段内测得的地理坐标系Y轴向的目标速度信息的滤波效果图；图10为应用EKF方法在30秒时间段内测得的地理坐标系Z轴向的目标速度信息的滤波效果图；图11为应用KF方法在30秒时间段内测得的地理坐标系Z辅向的目标速度信息的滤波效果图；图12为应用EMD-KF方法在30秒时间段内测得的地理坐标系Z轴向的目标速度信息的滤波效果图。

具体实施方式

具体实施方式一：本实施方式的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法，它的具体过程如下：

步骤一、获得当前采样时刻的量测数据序列，并利用卡尔曼预测方程，获得下一采样时刻的系统状态预测数据，并将所述下一采样时刻的系统状态预测数据和当前量测数据序列组合，获得组合后的数据序列；

步骤二、在步骤一获得的组合后的数据序列中的每相邻两个数据之间均匀分布多个插值点，生成滤波序列；

步骤三、利用EMD方法对步骤二生成的滤波序列进行分解，获得包含噪声的IMF，在滤波序列中剔除包含噪声的IMF，然后获得当前时刻的滤波值；

步骤四、将当前时刻的滤波值作为当前滤波结果显示；

步骤五、将当前时刻的滤波值作为系统当前状态的后验估计，在下一采样时刻来临时，结合下一采样时刻的量测数据，并利用卡尔曼方程组得到用于下一采样时刻的计算参数，然后，返回执行步骤一，将下一采样时刻作为当前采样时刻，实现机动目标的跟踪。

在系统运行过程中，若接收到来自外界的结束控制命令，则结束当前计算过程。

本发明综合运用了卡尔曼滤波(Kalman Filter-KF)理论和经验模态分解的方法，利用卡尔曼滤波方法对EMD滤波序列进行端点沿拓，并将EMD的滤波值代入卡尔曼递推方程中，有效地克服了EMD方法中“边界效应”以及针对强机动目标卡尔曼滤波的发散问题，具有较高的跟踪精度。本发明继承了EMD方法的自适应特性，因此对目标的机动方式、噪声的分布形式等条件并不敏感，具有较强的鲁棒性。

具体实施方式二：本实施方式是对实施方式一的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法的进一步说明，步骤一所述内容的过程如下：

在采样时刻t(n)，n＝1，2，…，滤波系统获取当前量测数据序列，且在各个采样时刻所获取的量测数据序列的长度均为N；

用{y(i) i＝1，2，…，N}表示t(n)时刻的量测数据序列，其中y(N)为t(n)时刻采样获得的系统状态量测数据，y(N-1)为t(n-1)时刻采样获得的系统状态量测数据，…，y(1)为t(n-N+1)时刻采样获得的系统状态量测数据；

所述卡尔曼预测方程为：

$\hat{x}(n+1|n)=F(n+1,n)\hat{x}\left(n\right|n),$

其中，n为当前采样时刻的序数，即当前采样时刻为t(n)；

F(n+1，n)为状态转移矩阵，它为已知量；

为用于预测的输入量，

为预测结果输出量；

在首次执行步骤一时，令

等于系统当前时刻的系统状态测量数据y(N)；在非首次执行步骤一时，令

等于上次执行步骤五时获得的所述用于下一采样时刻的计算参数；

通过卡尔曼预测方程，并根据

获得下一采样时刻t(n+1)的系统状态预测数据

将所述下一时刻的系统状态预测数据和原始量测数据序列组合，令{y(1)，y(2)，...，y(N)}表示原始量测数据序列，则组合后的数据序列为 $\{y\left(1\right),y\left(2\right),...,y\left(N\right),\hat{x}(n+1|n)\}.$

具体实施方式三：本实施方式是对实施方式一或二的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法的进一步说明，步骤二所述内容的具体过程为：

令表示步骤一获得的组合后的数据序列，其中y(1)、y(2)、……y(N)为步骤一获得的当前量测数据序列{y(i)，i＝1，2，…，N}中的N个系统状态量测数据，

为步骤一获得的下一采样时刻的系统状态预测数据；

在

中的每相邻两个数据之间分别均匀分布j个插值点，j为正整数，生成滤波序列{x(i)，i＝1，2，…，n×(j+1)+1}。

具体实施方式四：本实施方式是对实施方式一至三中的任意一种卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法的进一步说明，步骤三所述内容的具体过程为：

步骤三一、令{x(i)，i＝1，2，…，n×(j+1)+1}表示步骤二获得的滤波序列，利用EMD方法对该滤波序列进行分解，得到多个IMF，然后在所述多个IMF中获得含有噪声信号的IMF；

步骤三二、在滤波序列{x(i)，i＝1，2，…，n×(j+1)+1}中减去所述含有噪声信号的IMF，进而获得当前时刻的滤波值。

IMF为固有模态函数，满足以下两个条件的函数称为固有模态函数：

条件一：过零点和极值点个数相差不超过一；

条件二：上下包络线关于x轴对称。

其中，条件二不可能完全满足，只能近似满足。

具体实施方式五：本实施方式是对实施方式四的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法的进一步说明，步骤三一所述内容的具体过程为：

步骤三一一、令{x(i)，i＝1，2，…，n×(j+1)+1}表示步骤二获得的滤波序列，并将该滤波序列作为当前的待变换信号；初始化w＝0；

步骤三一二、令w＝w+1；获取当前的待变换信号的极大值点序列和极小值点序列，通过拟合获得所述极大值点序列的包络曲线以及所述极小值点序列的包络曲线；

步骤三一三、获得极大值点序列的包络曲线与极小值点序列的包络曲线的均值曲线m₁，然后开始进行迭代计算，并初始化h₀为待变换信号，所述迭代计算的公式为：

h_k＝h_k-1-m_k，k＝1，2，...

其中，h_k为第k次迭代计算的结果，h_k-1为第k-1次迭代计算的结果，k为迭代次数，m_k为h_k-1的极大值点序列的包络曲线与h_k-1的极小值点序列的包络曲线的均值曲线；

在所述迭代计算的过程中，当h_k满足IMF的定义条件时即停止迭代，并令c_w等于当前的迭代计算的结果，其中c_w为原始滤波序列的第w个IMF分量；

步骤三一四、计算r_w：

$r_w=\{x\left(i\right),i=\mathrm{1,2},\·\·\·,n\×(j+1)+1\}-\underset{q=1}{\overset{w}{\mathrm{\Σ}}}{c}_q;$

步骤三一五、判断r_w的极值点个数是否小于或等于2：若是，则执行步骤三一六；否则，返回执行步骤三一二；

步骤三一六、当前r_w即为信号残余；在w个IMF分量c₁、c₂....c_w中，获得含有噪声信号的IMF分量。

具体实施方式六：本实施方式是对实施方式五的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法的进一步说明，步骤三一六中所述的在w个IMF分量c₁、c₂....c_w中，获得含有噪声信号的IMF分量的具体过程为：

步骤A1、令 $\mathrm{Dif}\left(p\right)=\{x\left(i\right),i=\mathrm{1,2},\·\·\·,n\×(j+1)+1\}-\underset{q=1}{\overset{p}{\mathrm{\Σ}}}{c}_q,$ 并初始化p＝1；

步骤A2、计算当前Dif(p)，并判断当前Dif(p)所包含的极点个数是否小于预设的阈值：若是，则执行步骤A4；否则，执行步骤A3；

步骤A3、令p＝p+1，然后返回步骤A2；

步骤A4、在w个IMF分量c₁、c₂....c_w中，含有噪声信号的IMF分量为c₁、c₂....c_p。

需要指出的是，由于噪声分布和目标机动形式的不同，包含有噪声信号的IMF的阶次也是不同的。甚至对于噪声分布和机动形式相同的目标，由于采样数据长度的限制，包含有噪声信号的IMF的阶次也是在不断变化的，因此在滤波过程中需要引入决策机制加以判断。

具体实施方式七：本实施方式是对实施方式五的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法的进一步说明，步骤三一六中所述的在w个IMF分量c₁、c₂....c_w中，获得含有噪声信号的IMF分量的具体过程为：

分别依次判断w个IMF分量c₁、c₂....c_w的信号能量是否低于预设阈值，并将信号能量低于预设阈值的IMF分量作为含有噪声信号的IMF。

具体实施方式八：本实施方式是对实施方式四的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法的进一步说明，步骤三二所述内容的具体过程为：

在滤波序列{x(i)，i＝1，2，…，n×(j+1)+1}中减去所述含有噪声信号的IMF，获得一个新序列，则该新序列的第j+2个点的值即为当前时刻的滤波值。

具体实施方式九：本实施方式是对实施方式一至八中的任意一种卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法的进一步说明，步骤五所述内容的具体过程为：

令t(n)表示当前时刻，则将t(n)时刻的滤波值作为系统当前状态的后验估计代入如下的卡尔曼方程组中，即将步骤三中所得到的滤波值替换一下方程组中的

并结合t(n+1)时刻的测量值z(n+1)进行卡尔曼滤波估计：

$\hat{x}(n+1|n)=F(n+1,n)\hat{x}\left(n\right|n)$

$\hat{x}(n+1|n+1)=\hat{x}(n+1|n)+K(n+1)[z(n+1)-H(n+1)\hat{x}(n+1|n)]$

K(n+1)＝P(n+1|n)H^T(n+1)[H(n+1)P(n+1|n)H^T(n+1)+R(n+1)]^-1

P(n+1|n)＝F(n+1，n)P(n|n)F^T(n+1，n)+G(n+1，n)Q(n)G^T(n+1，n)

P(n+1|n+1)＝P(n+1|n)-K(n+1)H(n+1)P(n+1|n)

          ＝[I-K(n+1)H(n+1)]P(n+1|n)

即可获得用于下一采样时刻t(n+1)的计算参数该参数用于求取t(n+1)时刻的预测值

以下结合数字仿真算例阐述本发明的具体实施方式：

数字仿真是评判滤波跟踪算法性能的重要测试手段。本仿真环境及相关参数来自某雷达对运动目标的试验模拟数据。

运动目标作锥形螺旋机动。圆柱螺旋线的数学方程为：

$\left\{\begin{array}{c}{x}_t=\mathrm{\λt}\\ y_t=\mathrm{\ξ}\cos \mathrm{\ωt}\\ z_t=\mathrm{\ξ}\sin \mathrm{\ωt}\end{array}\right.$

其中，x_t、y_t和z_t分别为运动目标的空间三维坐标，t为时间变量，ξ为目标的螺旋机动半径，ω为其螺旋机动的角速度，λ为螺旋线轴线方向的速度。

由于目标的机动能力和大气密度、自身速度等因素有关，所以在不同高度的最大螺旋加速度a_max是不同的，数值由下式给出：

$a_{\max }=\left\{\begin{array}{cc}10g& H_t\≤10\mathrm{km}\\ 7g& H_t=19\mathrm{km}\\ 3g& H_t=23\mathrm{km}\\ g& H_t=27\mathrm{km}\end{array}\right.$

其中，H_t为目标所处高度，g为重力加速度。

本实施例中目标螺旋机动角速度ω＝3rad/s，初始螺旋线轴线方向的速度λ₀＝500m/s，ξ由a_MAX＝ω²ξ实时求出，弹道倾角在-50°-80°变化，雷达可以获取地理坐标系下目标的三维速度信息，且都含有均值为0，方差为10的随机噪声，采样周期0.1秒。仿真中将EMD-KF方法与当前被人们广为接受的卡尔曼方法及扩展卡尔曼方法的滤波效果做了对比，采用标准差作为评价指标，使得结果更加具有说服力。

在在采样时刻t(n)开始跟踪滤波，其中n＝1，2，…；将t(n)时刻的测量值和之前n-1个采样时刻的测量值按时间顺序排成数据序列{y(i)，i＝1，2，…，N}，其中y(N)为t(n)时刻采样获得的系统状态量测数据，y(N-1)为t(n-1)时刻采样获得的系统状态量测数据，…，y(1)为t(n-N+1)时刻采样获得的系统状态量测数据。根据测量值y(N)以及设定的状态转移矩阵F(n+1，n)，由卡尔曼理论的预测方程

可以推断出下一时刻系统的状态

与之前的量测数据构成新的数据序列此时数据序列长度为n+1，本实施例中，n＝5。

在新数据序列的每相邻两个数据之间均匀分布j个插值点，仿真中采用线性插值，j＝10，构成了可用于滤波的序列{x(i)，i＝1，2，…，n×(j+1)+1}。

应用EMD方法对序列{x(i)，i＝1，2，…，n×(j+1)+1}进行分解，其分解过程如下：

a)、通过下式迭代得到第一个IMF，公式中的x(t)是需要被变换的原信号，此处x(t)等于{x(i)，i＝1，2，…，n×(j+1)+1}，获取该信号的极大值点和极小值点的序列，用拟合的方法构造极值序列的包络曲线，m₁是它的上下包络线的均值，其余m_i是h_i-1上下包络线的均值。

$x\left(t\right)-m_1=h_1,$

$h_1-m_2=h_2,$

$h_{k-1}-m_k=h_k,$

$\⇒c_1=h_k$

h_k满足IMF的定义即可停止本次迭代，得到第一个IMF分量c₁。

b)、令r₁＝x(t)-c₁，对r₁重复上述迭代过程依次得到c₂，c₃…c_n，即：

$x\left(t\right)-c_1=r_1,$

$r_1-c_2=r_2,$

$r_{n-1}-c_n=r_n$

$\⇒x\left(t\right)=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{c}_j+r_n$

c)、当r_n的极值点小于或等于2时，则停止对x(t)的分解。

采用频率判断的方法，用原信号减去含有噪声的IMF后，得到最终的滤波信号，该信号序列中第j+2个点的滤波值即为当前时刻的滤波值。

将步骤三中所得到的滤波值作为系统当前状态的后验估计

并结合t(n+1)时刻的测量数据z(n+1)进行卡尔曼滤波估计：

$\hat{x}(n+1|n)=F(n+1,n)\hat{x}\left(n\right|n)$

$\hat{x}(n+1|n+1)=\hat{x}(n+1|n)+K(n+1)[z(n+1)-H(n+1)\hat{x}(n+1|n)]$

K(n+1)＝P(n+1|n)H^T(n+1)[H(n+1)P(n+1|n)H^T(n+1)+R(n+1)]^-1

P(n+1|n)＝F(n+1，n)P(n|n)F^T(n+1，n)+G(n+1，n)Q(n)G^T(n+1，n)

P(n+1|n+1)＝P(n+1|n)-K(n+1)H(n+1)P(n+1|n)

          ＝[I-K(n+1)H(n+1)]P(n+1|n)

经过迭代运算，可以得到参数

在t(n+1)采样时刻，用n+1代替n代入 $\hat{x}(n+1|n)=F(n+1,n)\hat{x}\left(n\right|n)$ 中，即通过 $\hat{x}(n+2|n+1)=F(n+2,n+1)\hat{x}(n+1|n+1)$ 获得t(n+2)采样时刻的系统状态预测数据

之后反复进行上述步骤，直到跟踪滤波过程的结束。

为了体现本发明所提设计方法的有效性和优越性，经过Monte Carlo随机仿真，采用标准差作为评价指标，将EMD-KF方法的滤波效果与当前被人们广为接受的卡尔曼方法(KF)及扩展卡尔曼方法(EKF)做了对比。表1为三种方法的滤波效果比较。

表1

	EMD-KF	KF	EKF
				X轴标准差	5.66	9.61	5.57
Y轴标准差	5.37	9.52	5.73
				Z轴标准差	5.51	8.84	5.47

应用EKF方法，在30秒时间段内，在地理坐标系X轴向、Y轴向和Z轴向获得的目标速度信息的滤波效果图分别为图4、图7和图10；

应用KF方法，在30秒时间段内，在地理坐标系X轴向、Y轴向和Z轴向获得的目标速度信息的滤波效果图分别为图5、图8和图11；

应用EMD-KF方法，在30秒时间段内，在地理坐标系X轴向、Y轴向和Z轴向获得的目标速度信息的滤波效果图分别为图6、图9和图12；

从表1及图4至图12中可以看出，卡尔曼滤波与经验模态分解结合的机动目标跟踪方法(EMD-KF)的滤波效果要优于KF方法，与EKF相当。而且由于EMD滤波的自适应特性，卡尔曼滤波与经验模态分解结合的滤波方法可以不需要事先获知噪声的统计特性，对噪声的分布、目标的机动形式也并不敏感，因此比EKF具有更好的适用性。另外，随着模型误差或目标机动能力的进一步加大，KF和EKF的滤波效果都出现了明显的降低甚至发散，而卡尔曼滤波与经验模态分解结合滤波方法的效果变化不大，体现了较强的鲁棒性。

Claims (8)

1.卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法，其特征在于它的具体过程如下：

步骤一、获得当前采样时刻的量测数据序列，并利用卡尔曼预测方程，获得下一采样时刻的系统状态预测数据，并将所述下一采样时刻的系统状态预测数据和当前量测数据序列组合，获得组合后的数据序列；

步骤二、在步骤一获得的组合后的数据序列中的每相邻两个数据之间均匀分布多个插值点，生成滤波序列；

步骤三、利用EMD方法对步骤二生成的滤波序列进行分解，获得包含噪声的IMF，在滤波序列中剔除包含噪声的IMF，然后获得当前时刻的滤波值；

步骤四、将当前时刻的滤波值作为当前滤波结果显示；

步骤五、将当前时刻的滤波值作为系统当前状态的后验估计，在下一采样时刻来临时，结合下一采样时刻的量测数据，并利用卡尔曼方程组得到用于下一采样时刻的计算参数，然后，返回执行步骤一，将下一采样时刻作为当前采样时刻，实现机动目标的跟踪；

步骤一所述内容的过程如下：

在采样时刻t(n)，n＝1，2，…，滤波系统获取当前量测数据序列，且在各个采样时刻所获取的量测数据序列的长度均为N；

用{y(i)，i＝1，2，…，N}表示t(n)时刻的量测数据序列，其中y(N)为t(n)时刻采样获得的系统状态量测数据，y(N-1)为t(n-1)时刻采样获得的系统状态量测数据，…，y(1)为t(n-N+1)时刻采样获得的系统状态量测数据；

所述卡尔曼预测方程为：

$\hat{x}(n+1|n)=F(n+1,n)\hat{x}\left(n\right|n),$

其中，n为当前采样时刻的序数，即当前采样时刻为t(n)；

F(n+1，n)为状态转移矩阵，它为已知量；

为用于预测的输入量，为预测结果输出量；

在首次执行步骤一时，令

等于系统当前时刻的系统状态测量数据y(N)；在非首次执行步骤一时，令等于上次执行步骤五时获得的所述用于下一采样时刻的计算参数；

通过卡尔曼预测方程，并根据

获得下一采样时刻t(n+1)的系统状态预测数据

将所述下一时刻的系统状态预测数据

和原始量测数据序列组合，令{y(1)，y(2)，...，y(N)}表示原始量测数据序列，则组合后的数据序列为 $\{y\left(1\right),y\left(2\right),...,y\left(N\right),\hat{x}(n+1|n)\}.$

2.根据权利要求1所述的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法，其特征在于步骤二所述内容的具体过程为：

令

表示步骤一获得的组合后的数据序列，其中y(1)、y(2)、……y(N)为步骤一获得的当前量测数据序列{y(i)，i＝1，2，…，N}中的N个系统状态量测数据，

为步骤一获得的下一采样时刻的系统状态预测数据；

在中的每相邻两个数据之间分别均匀分布j个插值点，j为正整数，生成滤波序列{x(i)，i＝1，2，…，n×(j+1)+1}。

3.根据权利要求2所述的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法，其特征在于步骤三所述内容的具体过程为：

步骤三一、令{x(i)，i＝1，2，…，n×(j+1)+1}表示步骤二获得的滤波序列，利用EMD方法对该滤波序列进行分解，得到多个IMF，然后在所述多个IMF中获得含有噪声信号的IMF；

步骤三二、在滤波序列{x(i)，i＝1，2，…，n×(j+1)+1}中减去所述含有噪声信号的IMF，进而获得当前时刻的滤波值。

4.根据权利要求3所述的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法，其特征在于步骤三一所述内容的具体过程为：

步骤三一一、令{x(i)，i＝1，2，…，n×(j+1)+1}表示步骤二获得的滤波序列，并将该滤波序列作为当前的待变换信号；初始化w＝0；

步骤三一二、令w＝w+1；获取当前的待变换信号的极大值点序列和极小值点序列，通过拟合获得所述极大值点序列的包络曲线以及所述极小值点序列的包络曲线；

步骤三一三、获得极大值点序列的包络曲线与极小值点序列的包络曲线的均值曲线m₁，然后开始进行迭代计算，并初始化h₀为待变换信号，所述迭代计算的公式为：

h_k＝h_k-1-m_k，k＝1，2，...

其中，h_k为第k次迭代计算的结果，h_k-1为第k-1次迭代计算的结果，k为迭代次数，m_k为h_k-1的极大值点序列的包络曲线与h_k-1的极小值点序列的包络曲线的均值曲线；

在所述迭代计算的过程中，当h_k满足IMF的定义条件时即停止迭代，并令c_w等于当前的迭代计算的结果，其中c_w为原始滤波序列的第w个IMF分量；

步骤三一四、计算r_w：

$r_w=\{x\left(i\right),i=\mathrm{1,2},\·\·\·,n\×(j+1)+1\}-\underset{q=1}{\overset{w}{\mathrm{\Σ}}}{c}_q;$

步骤三一五、判断r_w的极值点个数是否小于或等于2：若是，则执行步骤三一六；否则，返回执行步骤三一二；

步骤三一六、当前r_w即为信号残余；在w个IMF分量c₁、c₂....c_w中，获得含有噪声信号的IMF分量。

5.根据权利要求4所述的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法，其特征在于步骤三一六中所述的在w个IMF分量c₁、c₂...c_w中，获得含有噪声信号的IMF分量的具体过程为：

步骤A1、令 $\mathrm{Dif}\left(p\right)=\{x\left(i\right),i=\mathrm{1,2},\·\·\·,n\×(j+1)+1\}-\underset{q=1}{\overset{p}{\mathrm{\Σ}}}{c}_q,$ 并初始化p＝1；

步骤A2、计算当前Dif(p)，并判断当前Dif(p)所包含的极点个数是否小于预设的阈值：若是，则执行步骤A4；否则，执行步骤A3；

步骤A3、令p＝p+1，然后返回步骤A2；

步骤A4、在w个IMF分量c₁、c₂....c_w中，含有噪声信号的IMF分量为c₁、c₂....c_p。

6.根据权利要求4所述的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法，其特征在于步骤三一六中所述的在w个IMF分量c₁、c₂...c_w中，获得含有噪声信号的IMF分量的具体过程为：

分别依次判断w个IMF分量c₁、c₂....c_w的信号能量是否低于预设阈值，并将信号能量低于预设阈值的IMF分量作为含有噪声信号的IMF。

7.根据权利要求3所述的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法，其特征在于步骤三二所述内容的具体过程为：

在滤波序列{x(i)，i＝1，2，…，n×(j+1)+1}中减去所述含有噪声信号的IMF，获得一个新序列，则该新序列的第j+2个点的值即为当前时刻的滤波值。

8.根据权利要求1所述的卡尔曼滤波与经验模态分解有机结合的机动目标跟踪方法，其特征在于步骤五所述内容的具体过程为：

令t(n)表示当前时刻，则将t(n)时刻的滤波值作为系统当前状态的后验估计代入如下的卡尔曼方程组中，即将步骤三中所得到的滤波值替换一下方程组中的

并结合t(n+1)时刻的测量值z(n+1)进行卡尔曼滤波估计：

$\hat{x}(n+1|n)=F(n+1,n)\hat{x}\left(n\right|n)$

$\hat{x}(n+1|n+1)=\hat{x}(n+1|n)+K(n+1)[z(n+1)-H(n+1)\hat{x}(n+1|n)]$

K(n+1)＝P(n+1|n)H^T(n+1)[H(n+1)P(n+1|n)H^T(n+1)+R(n+1)]^-1

P(n+1|n)＝F(n+1，n)P(n|n)F^T(n+1，n)+G(n+1，n)Q(n)G^T(n+1，n)

P(n+1|n+1)＝P(n+1|n)-K(n+1)H(n+1)P(n+1|n)

          ＝[I-K(n+1)H(n+1)]P(n+1|n)

即获得用于下一采样时刻t(n+1)的计算参数

Images