基于学习的核偏最小二乘法优化扩展卡尔曼滤波

白晓波，邵景峰，和 征，田建刚

(1.西安工程大学管理学院，陕西 西安 710048; 2.陆军边海防学院信息与兵种教研室，陕西 西安 710108)

0 引 言

卡尔曼滤波(Kalman Filter, KF)于1960年由卡尔曼提出,为线性高斯问题提出了解决方案，至今被广泛应用于众多现实问题[1]。然而，很多科学或工程问题为非线性非高斯的，为了解决这类问题，扩展卡尔曼滤波(Extended Kalman Filter, EKF)[2]基于Taylor展开式，对非线性系统局部线性化，取得了较为理想的滤波效果，并得到广泛应用[3-5]。刘永进等[6]基于EKF对离散化后的故障模型设计扩展卡尔曼滤波器的残差产生器。文献[7]利用EKF估计水面和车辆的距离和速度。但是，EKF在针对强非线性系统时，滤波效果很不稳定，甚至发散。因此，众多学者对此进行了研究。文献[8]基于拉格朗日中值定理，以选择更好的线性展开点和雅可比矩阵的计算点，对EKF进行改进，改进算法的滤波精度有所提高，但是却较大地增加了运算量，进而影响算法实时性。刘毛毛等[9]基于多新息理论对EKF的估计值进行修正，在算法运算量增加不大的情况下，提高了EKF的滤波精度，但是新息数的选择受限于现实系统，而新息数的选择影响滤波精度、算法收敛性和实时性。文献[10]通过计算加权高斯积分的基本容积点及其对应权值对EKF进行改进。刘国海等[11]提出基于插值公式的EKF改进算法。宁倩慧等[12]在考虑斜距误差和多普勒误差相关性的情况下将多普勒量测包含于扩展卡尔曼滤波算法，从而提高了目标跟踪的精度。吴汉洲等[13]提出基于多项式拟合的EKF改进算法，以提高滤波精度，并在一定程度上提高了算法实时性。其他如文献[14-15]利用迭代的方法提高滤波精度，但是却较大地增加了算法运算量，降低了算法实时性。

国内外学者都提出了扩展卡尔曼滤波的改进算法，但是滤波效果与系统方程、过程噪声方差和量测噪声方差紧密相关，且这些参数具有时变特性，使得算法的应用有其局限性。因此，研究系统参数不确定、过程方差和量测方差信息不准确导致的EKF估计随时序误差较大的问题，具有较强的理论意义和工程实践意义。为此，本文提出基于核偏最小二乘法优化EKF的滤波方法KPLS-EKF。

1 扩展卡尔曼滤波与问题描述

扩展卡尔曼滤波EKF的核心思想是，利用非线性函数泰勒级数展开线性化方法近似表示原系统的状态方程和量测方程。

1.1 扩展卡尔曼滤波

扩展卡尔曼滤波的状态方程和量测方程表示如下：

(1)

式中，xk为系统状态矩阵，zk为量测矩阵，uk为系统输入量，f(·)为非线性映射函数，利用上一时刻k-1时的估计状态计算当前时刻k时的预测状态；h(·)以预测的状态计算预测的测量值；wk和vk分别为系统状态噪声矩阵和量测噪声矩阵。wk和vk相互独立，均值和协方差如下：

(2)

式中，δ(k-j)为Kronecker函数。EKF算法如下所示[16]。

预测方程：

(3)

(4)

更新方程：

(5)

(6)

(7)

(8)

Pk|k=(I-KkHk)Pk|k-1

(9)

式中：

1.2 问题描述

EKF随时序变化的滤波效果出现较大偏差。这是由于EKF在非线性函数泰勒展开时仅用了一阶项，而忽略了高阶项，同时在EKF中构建的系统状态方程和量测方程的参数的不确定和噪声信息的不准确，即模型参数和噪声随时序变化，为滤波效果带来了较大误差。如图1所示。

图1 扩展卡尔曼滤波估计误差

在图1中表示了误差随时序变化的趋势。因此，为了解决由于系统参数不确定和噪声信息不准确引起的估计偏差随时序增大的问题，在EKF中融入KPLS学习算法，以抑制其误差增大。

2 基于KPLS优化EKF

2.1 核偏最小二乘法KPLS

为了解决非线性系统的预测问题，Trejo和Rosipal教授[17]将核函数引入PLS，提出核偏最小二乘法(Kernel Partial Least Square, KPLS)，其基本原理是，通过非线性核函数将输入空间映射到高维特征空间，在特征空间中构建线性偏最小二乘回归，从而实现原始输入空间的非线性建模。核偏最小二乘算法的详细步骤如下[18]。

步骤1计算核矩阵：

(10)

利用高斯核函数计算矩阵元素，高斯核函数表示如下：

(11)

式中σl(l=1,2,…,p)表示高斯核函数宽度。

步骤2对核矩阵K中心化处理：

(12)

步骤3随机初始化Y的得分向量u。

步骤4计算特征空间中X的得分向量th并进行归一化：

(13)

步骤5计算Yh的权值向量ch：

(14)

步骤6计算Yh的得分向量uh并归一化：

(15)

步骤7重复步骤4～步骤6，直到th收敛。

步骤8将矩阵K、 Y缩小，重复步骤3～步骤7，取得p个t、 u：

(16)

(17)

以下为训练样本的拟合公式：

(18)

预测样本拟合公式为：

(19)

Kt=Φ(xnew)Φ(x)T

(20)

2.2 KPLS-EKF算法

2.2.1 KPLS优化EKF理论分析

根据式(3)、式(5)和式(8)可知，在标准的EKF滤波算法中，k时的预测，只利用了线性化后k-1时刻的状态估计，而忽略了k-2，k-3，…，k-p时的状态估计信息，这样隐含的系统参数时变信息就不能在以后的状态估计中体现。

KPLS具有较强的泛化性能和抗噪能力，且在预测精度和训练时间与SVM(Support Vector Machine)相当。而KPLS基于历史采样数据，实现从原始输入空间到输出的非线性建模和预测，且不用关注非线性系统的状态方程和过程方差等信息。因此，可以利用KPLS建立量测与EKF收敛估计之间的核模型，再以此模型计算k时的预测值，并和EKF预测值合成，将此合成值作为k时状态估计的重要参数。

基于以上分析，根据EKF对非线性系统状态的估计过程，基于KPLS对动态非线性系统的模型重新定义如下：

Y=X·B

(21)

(22)

式中Y为系统状态各分量在k时的状态预测，X为系统各分量xl在k-1,k-2,k-3,…,k-p时刻的量测矩阵，B为模型的参数向量，由核偏最小二乘法辨识得到。

2.2.2 基于滑动窗口更新KPLS预测模型

基本思想是设定固定长度的样本作为训练样本，并使窗口随着新样本的采集不断向前滑动，以替换旧样本，从而模型能够定期更新，实现KPLS实时预测。具体过程如下。

步骤1对于变量X∈Rn×p，因变量Y∈Rn，自变量个数为p, n表示样本个数。其输入构成一个n×p的矩阵，用公式(6)和公式(7)处理后，得核矩阵K。因此，设定窗口长度为n，滑动距离为q。

图2 滑动窗口更新采样数据过程

2.2.3 高斯核函数宽度参数σ的取值

在KPLS中，重点在于核函数中的宽度参数σ的取值，其值影响KPLS的性能。σ较大时，主要应用于平缓变化的样本，具有更优的泛化能力；σ较小时，主要应用于剧烈变化的样本。在多输入变量的情况下，各变量的宽度参数σ不同，且随时序变化。通常σ的取值可参照小波变换中尺度变化规律，多尺度小波核有近似正交的特点，适用于连续函数的拟合及数字信号的局部分析、信噪分离和突变信号的检测。具体方法如下。

σi=2iσ, i=0,1,2,…

(23)

2.2.4 KPLS-EKF量测噪声协方差估计

为了增强KPLS-EKF的适应性，需要解决噪声信息不准确的问题。通常，在标准EKF中，滤波增益与系统状态估计协方差和量测噪声协方差都有关，而量测噪声对滤波效果有更大的影响。因此，基于Sage-Husa算法对量测噪声协方差进行估计，而忽略系统状态估计协方差的估计，以提高算法效率。量测噪声协方差的估计如下所示[20]。

R(k)= dk{[1-H(k)K(k-1)]V(k)VT(k)·

[1-H(k)K(k-1)]T+H(k)P(k-1)HT(k)}+

(1-dk)R(k-1)

(24)

式中，dk为遗忘因子：

(25)

在滤波过程中，量测噪声协方差并非需要及时更新。那么就需要明确量测噪声协方差的更新条件。通常，采用滤波收敛性作为判据，判据表示方法如下：

V(k)VT(k)H(k)P(k|k-1)HT(k)+R(k)

(26)

若式(26)为真，则滤波收敛，无需更新量测协方差，令R(k)=R(k-1)；否则，利用式(24)更新量测协方差。

2.2.5 KPLS-EKF算法

基于2.2.1节～2.2.4节的内容，设计KPLS-EKF算法流程如图3所示。

KPLS-EKF的详细步骤如下。

步骤1采样与KPLS建模阶段。

(a)利用标准EKF对系统状态进行初始估计，以形成历史数据。

(b)用式(26)判定EKF是否收敛，若收敛，则执行步骤(c)；否则，返回步骤(a)进行下一时刻的预测。

(c)将量测数据与状态估计值采样。

(d)判断样本数是否等于N，若等于N，则执行步骤(e)；否则，返回步骤(a)进行下一时刻系统状态估计。

(e)以量测数据作为输入矩阵，状态估计值矩阵作为输出，利用式(10)～式(19)初始化KPLS模型。

图3 KPLS-EKF算法流程图

步骤2KPLS预测与模型更新。

(g)KPLS窗口向前滑动，若滑动距离d=q，执行步骤(h)；否则，KPLS模型不变，返回步骤(f)。

(h)生成随机数i, i∈[1,5]，利用公式(23)取核参数σ，再更新KPLS模型，返回步骤(f)。

步骤3EKF预测与更新。利用EKF的公式(3)预测得k时状态xk|k-1。

(27)

式中ρKPLS表示KPLS预测信度，ρEKF表示EKF预测信度，理想状态下ρKPLS=ρEKF=1/2，即KPLS和EKF预测的信度相同。

步骤5利用EKF的公式(5)～公式(8)进行状态估计得最终系统状态xk|k，并将其作为KPLS模型更新的样本。

步骤6利用公式(9)更新协方差矩阵Pk|k。

步骤7若式(26)成立，则系统量测噪声协方差不变，执行R(k)=R(k-1)；否则，用式(24)更新R(k)。

在以上KPLS-EKF算法中，融入了基于滑动窗口的KPLS，通过时序采样、训练与学习，构建量测与状态的收敛估计的核矩阵，而无需建立具体的系统状态模型，可以避免由于系统模型参数的不准确导致的估计误差。另外，通过滤波收敛性判定量测协方差的更新，提高KPLS-EKF算法在量测噪声信息不准确的情况下的滤波精度。

2.3 KPLS-EKF时间复杂度

从理论上分析，KPLS优化EKF的可行性KPLS-EKF是在EKF的基础上融入KPLS，而EKF本身时间复杂度较低，具有较高的实时性。在融入基于学习的KPLS后，整个算法的时间复杂度必定高于标准EKF的时间复杂度，因此针对KPLS-EKF算法的详细过程进行分析，有利于选取合适的参数，在取得预期滤波精度的情况下，实时性也能够得到保证。由图3及其详细步骤，可以将KPLS-EKF算法分为如下3部分。

1)数据采集和KPLS模型初始化。该阶段作为模型的学习和训练阶段，其运算量取决于2方面，一是采样的样本数N，二是非线性迭代求解T和U时，设置的终止迭代次数t和收敛判据，以及输入空间自变量个数m，若不能满足收敛判据，达到迭代终止条件后则终止。那么，最坏情况下，学习训练KPLS模型的时间复杂度为O(N2+m×t)。但是，该阶段作为模型学习与训练的必要时间开销，可以不作为系统状态估计时的时间开销。

3)KPLS模型更新阶段。该阶段从本质上分析，是在达到模型更新条件时对核矩阵进行更新，整个过程和核矩阵初始建模没有任何区别，故其时间复杂度仍然表示为O(N2+m×t)。

综合以上分析，在系统状态实时估计阶段，KPLS-EKF的时间复杂度表示为O(2N2+M3+m×t)。在实际应用时，t、N和M取值都较小，一般在10以内，而即使迭代次数t=100，最坏情况下其迭代运算量仅为1000。因此，综合分析，基于学习的KPLS优化EKF解决非线性系统的状态估计，在理论上具有可行性。

3 实验仿真与分析

利用实验仿真的方法重点探究以下3个问题。

1)量测与估计的样本数N的取值。N的大小影响到2个方面。一是算法时间复杂度。N越大，算法运算量越大；反之，运算量越小。二是合适的样本数量有利于构建准确的KPLS核矩阵，从而提高算法的收敛性和预测准确度。但是运算效率与收敛性和精度之间是个矛盾的问题，因此，这里需要找到合适的样本数N，在效率与性能之间找到平衡。

2)KPLS-EKF算法收敛性分析。主要从实验仿真的角度分析KPLS-EKF在经过多少次迭代以后，估计值与真值的误差趋向于稳定，或者相较于固定值合理波动。

3)KPLS-EKF性能分析。重点和其它算法相比，其滤波精度如何，主要通过均方根误差(Root Mean Square Error, RMSE)作对比分析。由于每次实验具有随机性，因此，每次取均方根误差的均值作为对比依据。

以上实验的运行环境为Windows XP系统，Intel(R) Core(TM) i3-2130 CPU @ 3.40 GHz & 3.39 GHz, 4 GB内存，500 GB硬盘，仿真软件为Matlab R2012a。实验的非线性系统状态方程如下：

Xk=f(xk-1,uk,wk)=F×Xk-1+G×sqrt(Q)×rand(2,1)

式中：

F为状态转移矩阵，G为过程噪声驱动矩阵，Q为过程噪声方差；T为雷达扫描周期，观测站的位置为Xstation=[x0,y0]。量测方程如下：

Z(k)=Dist(X(:,k),Xstation)+sqrt(R)·rand

式中Dist(·)表示目标与观测站的距离，R为量测方差。

系统的初始化数据为T=1，总采样次数N=200/T,初始状态为X(:,1)=[-100,10,200,20]，-100为x方向位置，10为x方向速度，200为y方向位置，20为y方向速度。观测站的位置x0=200, y0=300，即Xstation=[200,300]。

3.1 KPLS-EKF样本数量选择

通过多次实验仿真的方法，分析样本数量的选择对滤波精度的影响。具体方法如下。

分别取样本数N=2,3,…,20建立KPSL核矩阵，并分别进行50次仿真实验。在此基础上记录下选取每个样本数时滤波的RMSE，并计算50次实验的RMSE均值，计算方法如下：

(28)

式中M为取不同的样本数时实验仿真次数，N为建立KPLS核矩阵时的样本数。下面所有的RMSE均值计算都基于公式(28)。以X方向的速度估计为例，实验结果如图4所示。

(a) N=2

(b) N=3

(c) N=4

(d) N=5

在图4中，滤波的起始为KPLS-EKF采样和学习阶段，在达到训练样本以后，有相应的滤波结果输出。在后续的滤波阶段，KPLS-EKF的值更接近于EKF，即训练样本为2时，能够较快地建立KPLS核矩阵，但是其滤波效果更接近EKF的估计值，而在N=4时，其效果更加显著。而到了N=5时，效果并不明显，这是因为过多历史数据的采样反而会抑制基于训练样本构建的核模型带来的有效增益。因此，选取合适的训练样本数显得尤为重要。

基于以上4项实验结果，初步得出如下结论。即样本数的增加，KPLS-EKF的滤波结果逐步靠近真值，虽然较EKF能够较快收敛于真值附近，却需要样本收集和训练为代价。为了进一步分析KPLS-EKF中训练样本对滤波结果的影响，考虑实验的随机性，取样本数N=2,3,…,20，对每个取值进行50次实验，并利用式(28)计算在X方向上速度估计值和真值的RMSE均值，得出如图5所示的结果。

训练样本数N∈[3,8] 之间为宜。这是因为在图5中，整体上随着样本数N 的增加，RMSE均值减小，滤波精度逐步提高。但是，随着样本数的增加，RMSE均值稳定在一定范围内，呈一定的波动性，主要是因为实验时量测与估计值都具有随机性，其波动值在合理范围内，且过多地增加训练样本并不能有效地增加算法的滤波精度。另外，根据2.3节的分析，过多的训练样本必然带来时间复杂度的增加，影响算法的实时性。

图5 不同样本数下KPLS-EKF的RMSE均值

3.2 KPLS-EKF收敛性分析

为了分析KPLS-EKF的收敛性，根据3.1节的实验结论，设定样本数量N=4，观察系统状态2个参数x方向速度和y方向速度的收敛情况。实验结果如图6所示。

(a) X方向速度

(b) Y方向速度

图6所示的实验结果与图4中样本数N=4时的结果接近。在KPLS-EKF样本的搜集阶段，以EKF的估计值作为滤波结果，在达到样本数以后，能够较快地收敛于真值附近，用了近20步，EKF为40步。以上实验结果具有一定的随机性，因此，进行了50次实验求KPLS-EKF、标准EKF和改进MI-EKF[21]，在x和y的2个方向速度的收敛步数均值。结果如表1所示。

表1中的p表示新息数，N 为训练样本数。

从表1中的2个参数Vx、 Vy在不同算法的不同参数下的收敛步数可以得出以下结论：

1)KPLS-EKF算法收敛速度和改进MI-EKF算法的接近，均优于标准EKF算法。但在整体上，KPLS-EKF的收敛速度略好于改进MI-EKF算法。

2)通过增加样本数，KPLS-EKF算法能够提高收敛速度，但是，增加到一定程度后，其收敛性并不能显著提高，这是因为合适的采样数据能够对滤波产生积极的优化效果，然而，过多的训练数据对最终估计值产生了累积干扰。所以，需要在有效优化和抑制累积之间找到相对平衡点。

表1 3种滤波算法同一参数收敛步数均值

序号滤波算法参数收敛步数VxVy1EKF42482KPLS-EKF(N=2)40453KPLS-EKF(N=3)35374KPLS-EKF(N=4)25235KPLS-EKF(N=5)24236改进MI-EKF(p=2)38407改进MI-EKF(p=3)35368改进MI-EKF(p=4)30339改进MI-EKF(p=5)2930

3.3 KPLS-EKF性能分析

为了进一步说明KPLS-EKF算法的有效性，本文选取训练样本数为4，将KPLS-EKF算法与标准EKF算法、MI-EKF算法[9]和文献[21]的改进MI-EKF算法进行RMSE均值对比，以分析其滤波精度和运行效率。

(a) X方向速度

(b) Y方向速度

根据图7中4种滤波算法的对比结果，初步得出结论为：KPLS-EKF在滤波的初始阶段，未采集到足够的收敛训练样本，使用EKF的估计值，故其结果和EKF相同；在达到训练样本数以后，KPLS-EKF的估计值相比其它几种滤波算法，更靠近真实值。但是，考虑实验仿真过程存在随机性，重复进行50次实验，设过程噪声方差Q=5×diag([0.5,1])，观测噪声方差R=5+rand([1,1])，利用式(28)计算4种算法的估计值和真值的均方根误差(RMSE)均值。实验结果如表2所示。

表2 4种算法滤波性能对比

序号滤波算法RMSE均值VxVy运行时间均值/s1KPLS-EKF(N=2)8.652611.76180.02252KPLS-EKF(N=3)4.24327.35240.02313KPLS-EKF(N=4)2.16844.27760.02384EKF9.109112.21830.02125MI-EKF(p=2)7.368210.47740.02196MI-EKF(p=3)5.56268.67180.02267MI-EKF(p=4)2.38515.49430.02358改进MI-EKF(p=2)7.175910.28510.02219改进MI-EKF(p=3)5.36258.47170.022410改进MI-EKF(p=4)2.27444.38860.0236

表2中参数p、 N 的含义同表1。

由表2的实验数据得出结论如下：

1)KPLS-EKF在算法效率上较MI-EKF和改进的MI-EKF略有降低，是因为在滤波过程中，需要不断更新核矩阵，增加了算法运算量，但是增加量并不明显，当样本数2N4 时，运行时间增加值在0.0001 s～0.0007 s之间，增加幅度在可接受范围内。

2)KPLS-EKF在训练样本数N>2时，RMSE均值均小于EKF、 MI-EKF和改进的MI-EKF。在样本数N=3时，KPLS-EKF在X方向速度的RMSE均值为4.2432, Y方向速度的RMSE均值为7.3524，效果优于其它3种算法。而到了样本数为4，MI-EKF和改进的MI-EKF的新息数为4时，这3种算法的滤波效果相近。

4 结束语

扩展卡尔曼滤波在对非线性系统状态估计时，需要状态方程准确的参数、量测与噪声方差信息，且这些参数信息存在时序变化的问题，所以通常建立的系统状态方程，难以适应时变系统的状态估计。因此，将KPLS融入扩展卡尔曼滤波，即通过学习算法，建立量测与EKF收敛估计值之间的预测模型，从而避免了系统状态方程、量测与过程方差信息不准确导致EKF估计值随时序发散的问题。然而，在算法的初始采样阶段，其滤波精度并未有效改善，另外，在固定的训练样本数下导致算法运行效率低于MI-EKF和改进的MI-EKF，下一步应该在此基础上研究提高算法效率的方法。