一种基于多胞型微分包含的非线性滤波方法

摘要

本发明涉及一种基于多胞型微分包含的非线性滤波方法，属于系统滤波与控制技术领域。本方法将非线性滤波误差系统采用PLDI模型来描述，从而将非线性滤波算法设计问题转换为线性不确定系统鲁棒滤波算法设计问题；再利用混合鲁棒H2/H∞滤波方法，设计估计误差校正量求解的动态方程；然后结合EKF一步预测方程，设计离散型非线性滤波方程；将其用于非线性离散系统中，实时获取其状态估计。本发明简化了非线性滤波设计，不需要实时更新滤波增益，且在实现过程中不需要实时计算雅可比矩阵，极大地降低了计算量，有效地提高了非线性滤波的实时性；适用于非线性滤波器设计。

说明书

技术领域

本发明涉及一种基于多胞型微分包含的非线性滤波方法，属于系统滤波与 控制技术领域。

背景技术

对于非线性滤波问题，在理论上很难找到严格的最优解，一般采用近似方 法来求解。根据对系统非线性处理方法的不同，非线性滤波方法主要可分为三 大类，第一类是函数近似的方法，采用泰勒级数展开或者插值多项式展开的方 法对非线性函数进行近似；第二类是确定性采样方法，即对非线性概率密度函 数进行近似的方法；第三类是基于蒙特卡洛仿真的方法。其中典型代表算法且 应用最广泛的是扩展卡尔曼滤波方法(EKF)。

EKF采用线性化近似的方法，对非线性函数进行泰特级数展开，并保留至 一阶，从而将非线性系统转换为线性系统，进而采用卡尔曼滤波方法实现对非 线性系统的状态估计。EKF为迭代算法，虽然简单易于实现，但由于在线性化 过程中引入了模型误差，使得滤波精度下降，甚至会出现滤波发散现象，且在 滤波过程中需要实时求解雅克比矩阵，计算复杂，尤其在高维非线性滤波求解 中，很容易出现求解困难或错误等问题。针对EKF存在的不足，学者们提出了 很多改进方法，如高阶截断EKF，变增益EKF、加权EKF、基于神经网络的EKF 方法。改进的EKF方法虽然在一定程度上提高了系统滤波的精度，但同时增加 了滤波算法的计算量，且在实质上仍是采用线性化近似的方法来实现对非线性 状态的估计，本质上并未克服EKF算法中的不足。

一般情况下，近似非线性的概率密度函数比近似非线性函数要容易得多， 基于该思想产生了近似求解非线性密度函数的非线性滤波方法，典型的代表算 法为Unscented卡尔曼滤波(UKF)和粒子滤波方法。UKF的核心是UT变换，与 EKF相比，UKF具有很好的一些特性：(1)利用UT变换对非线性函数的概率密 度分布进行近似，避免了线性化过程中引入的模型误差；(2)非线性分布统计量 的计算精度至少达到二阶；(3)不需要求解雅克比矩阵，在计算量上也没有太大 的增加。但综合考虑滤波性能与计算量时，UKF与EKF相比较，EKF是较好的 选择。粒子滤波是一种基于蒙特卡洛方法实现Bayes估计的一种滤波方法，理 论上，粒子滤波是最优的，且适用于非高斯噪声情况，但存在粒子退化等主要 问题。针对粒子退化问题，学者提出了很多改进算法，如辅助粒子滤波、正则 粒子滤波、高斯粒子滤波、Unscented粒子滤波算法、Rao-Blackwellized粒子滤 波等，在一定程度上改善了粒子退化现象，但仍存在计算量大、实时性差、滤 波稳定性差等问题。

微分包含理论通过采用全局线性化方法，将非线性系统采用线性微分包含 (LDI)模型来描述，原非线性系统为LDI系统(LDIs)的一个子集，虽具有一定的 保守性，但由于LDIs的线性特性，为简化非线性系统控制及滤波算法的设计提 供了新的思路。

发明内容

本发明的目的是为解决一般非线性系统的滤波算法复杂、计算量大的问题， 提出了一种基于多胞型微分包含技术的非线性滤波方法，该方法将非线性滤波 误差系统采用PLDI模型来描述，从而将非线性滤波算法设计问题转换为线性不 确定系统鲁棒滤波算法设计问题。

本发明的一种基于PLDI的非线性滤波方法，具体包括以下步骤：

步骤1，建立非线性滤波误差系统的PLDI描述模型。

非线性离散系统表达式为：

x_k＝f(x_k-1)+w_1k-1+v_1k-1                       (1)

y_k-1＝h(x_k-1)+w_2k-1+v_2k-1

其中，x_k-1∈Rⁿ为系统k-1时刻的状态，y_k-1∈R^m为k-1时刻的系统测量输出，n、m 分别表示系统状态变量的维数和测量输出的维数；f(x_k-1)、h(x_k-1)为非线性连续 可微函数；w_1k-1、v_1k-1为k-1时刻的系统过程噪声，为零均值不相关的白噪声；w_2k-1、 v_2k-1为k-1时刻的测量噪声，为能量有限的不相关噪声。

非线性离散系统在k时刻的状态估计值为定义状态估计误差和测量估计 误差分别为：

$\Delta x_k=x_k-{\hat{x}}_k,$${\mathrm{\Delta y}}_{k-1}=y_{k-1}-{\hat{y}}_{k-1}$

得到滤波误差系统模型为：

$\Delta x_k=f\left(x_{k-1}\right)-f\left({\hat{x}}_{k-1}\right)+B{w}_{k-1}+{\mathrm{Bv}}_{k-1}---\left(2\right)$

$\Delta y_{k-1}=h\left(x_{k-1}\right)-h\left({\hat{x}}_{k-1}\right)+{\mathrm{Dw}}_{k-1}+{\mathrm{Dv}}_{k-1}$

其中，B=[I_n×n 0_n×m]，D=[0_m×n I_m×m]，$w_{k-1}=\left(\begin{array}{c}{w}_{1k-1}\\ w_{2k-1}\end{array}\right),$$v_{k-1}=\left(\begin{array}{c}{v}_{1k-1}\\ v_{2k-1}\end{array}\right);$

若x_k∈Ω₁，y_k-1∈Ω₂，Ω₁、Ω₂为紧集，则存在PLDI描述模型：

Δx_k＝AΔx_k-1+Bw_k-1+Bv_k-1             (3)

Δy_k-1＝CΔx_k-1+Dw_k-1+Dv_k-1

包含上述滤波误差系统模型(2)。其中，

$(A,C)\in \left\{(A,C)\right|(A,C)=\underset{j=1}{\overset{l}{\Sigma }}{\lambda }_j(A_j,C_j),0\le {\lambda }_j\le 1,\underset{j=1}{\overset{l}{\Sigma }}{\lambda }_j=1\}$

(A_j,C_j)PLDI描述模型的第j个顶点，l为PLDI的顶点数。

步骤2，利用混合鲁棒H2/H∞滤波方法，设计估计误差校正量求解的动态 方程。

由于PLDI描述模型(3)包含原滤波误差系统模型(2)，根据对父系统适用的 滤波器对子系统同样适用的特性，结合离散型鲁棒H2/H∞滤波方法，得到非线 性离散系统的估计误差误差校正量的动态方程为：

$\Delta {\hat{x}}_{\mathrm{mk}}=A_F\Delta {\hat{x}}_{\mathrm{mk}-1}+B_F\Delta y_{k-1}$(4)

$\Delta {\hat{x}}_{k-1}=C_F\Delta {\hat{x}}_{\mathrm{mk}-1}$

其中，为k-1时刻的估计误差校正量，为k-1时刻的估计误差校正量的 预测，A_F、B_F、C_F为滤波系数：

A_F＝G^-1S_A、B_F=G^-1S_B、C_F=S_C                 (5)

其中，矩阵S_A、S_B、S_C、G通过求解下述优化问题得到：给定常数γ＞0，对于 正定矩阵变量P_11j、P_22j、Q_11j、Q_22j和矩阵变量P_12j、Q_12j、R、S、G、S_A、S_B、 S_C，寻求矩阵Z的最优上界，使得不等式(i)、(ii)、(iii)成立。

$\left(i\right)\left(\begin{array}{ccc}{Q}_{11j}-R^T-R& Q_{12j}-G-S& R^{T}B+S_{B}D\\ *& Q_{22j}-G-G^T& S^{T}B+S_{B}D\\ *& *& -Z\end{array}\right)<0$

$\left(\mathrm{ii}\right)\left(\begin{array}{ccccc}{Q}_{11j}-R^T-R& Q_{12j}-G-S& R^T{A}_j+S_B{C}_j& S_A& 0_{n\times n}\\ *& Q_{22j}-G-G^T& S^T{A}_j+S_B{C}_j& S_A& 0_{n\times n}\\ *& *& -Q_{11j}& -Q_{12j}& I_{n\times n}\\ *& *& *& -Q_{22j}& -S_C\\ *& *& *& *& -I_{n\times n}\end{array}\right)<0$

$\left(\mathrm{iii}\right)\left(\begin{array}{cccccc}{P}_{11j}-R^T-R& P_{12j}-G-S& R^T{A}_j+S_B{C}_j& S_A& R^{T}B+S_{B}D& 0_{n\times n}\\ *& P_{22j}-G-G^T& S^T{A}_j+S_B{C}_j& S_A& S^{T}B+S_{B}D& 0_{n\times n}\\ *& *& -P_{11j}& -P_{12j}& 0_{n\times (n+m)}& I_{n\times n}\\ *& *& *& -P_{22j}& 0_{n\times (n+m)}& -S_C\\ *& *& *& *& -\gamma I_{(n+m)\times (n+m)}& 0_{(n+m)\times n}\\ *& *& *& *& *& -I_{n\times n}\end{array}\right)<0$

其中，j=1,2,…,l,滤波估计误差协方差矩阵满足tr(·) 表示矩阵的迹。

步骤3，结合EKF一步预测方程，设计离散型非线性滤波方程。

结合EKF的一步预测方程及上述滤波误差校正量求解的动态方程(4)，得到 离散型非线性滤波方程为：

状态一步预测：${\hat{x}}_{k/k-1}=f\left({\hat{x}}_{k-1}\right)---\left(6\right)$

状态估计误差校正量预测：$\Delta {\hat{x}}_{\mathrm{mk}}=A_F\Delta {\hat{x}}_{\mathrm{mk}-1}+B_F[y_{k-1}-h\left({\hat{x}}_{k-1/k-2}\right)]---\left(7\right)$

状态估计误差校正量：$\Delta {\hat{x}}_k=C_F\Delta {\hat{x}}_{\mathrm{mk}}---\left(8\right)$

系统的状态估计${\hat{x}}_k={\hat{x}}_{k/k-1}+\Delta {\hat{x}}_k---\left(9\right)$

步骤4，将步骤3得到的离散型非线性滤波方程用于步骤1所述的非线性离 散系统中，实时获取其状态估计。

具体实现过程为：

步骤4.1，给定初始状态估计值及误差校正量预测的初始值

步骤4.2，对于k＞0时刻，根据k-1时刻的系统状态估计结合状态一步 预测方程获取k时刻的系统状态一步预测

步骤4.3，利用过程中得到的k-1时刻系统观测值y_k-1、状态一步预测及 状态估计误差校正量预测结合状态估计误差校正量预测方程，得到的k 时刻状态估计误差校正量预测其中k=1时刻，状态估计误差校正量预测 用步骤4.1中给定的初始值；

步骤4.4，结合状态估计误差校正量方程，获取k时刻状态估计误差校正量 并用对步骤4.2得到的k时刻系统状态一步预测值进行校正，最终 获取k时刻的系统状态估计

有益效果

本发明方法重点考虑非线性滤波器算法设计问题，基于PLDI技术，将非线 性滤波误差系统采用PLDIs模型来描述，从而将非线性滤波器设计问题转换为 线性不确定系统鲁棒滤波器设计问题，简化了非线性滤波的设计方法；在此基 础上，利用混合鲁棒H2/H∞滤波技术，获取求解滤波误差校正量的动态方程， 并对EKF算法得到的系统状态估计的一步预测值进行校正，由于求解滤波误差 校正量的动态方程中的滤波系数为常量，不需要实时更新滤波增益，且在实现 过程中不需要实时计算雅可比矩阵，极大地降低了算法的计算量，有效地提高 了非线性滤波算法的实时性；适用于非线性滤波器设计。

附图说明

图1为本发明中的基于多胞型微分包含新的非线性滤波方法原理示意图；

图2为具体实施方式中高度绝对估计误差对比曲线图；

图3为具体实施方式中速度绝对估计误差对比曲线图；

图4为具体实施方式中弹道系数绝对估计误差对比曲线图。

具体实施方式

为了更好地说明本发明的目的和优点，下面结合附图和实施例对本发明作 进一步说明。

一般情况下，非线性滤波算法复杂，计算量大，难以满足非线性状态估计 实时性的要求，非线性滤波算法的复杂性与计算量是非线性滤波算法实时性的 重要指标，对算法能否在实际工程应用中得以实现有着重要的影响作用。为比 较本发明中提出的非线性滤波算法和EKF的实时性，以滤波算法的计算量的大 小表征算法的复杂性和实时性指标，得到本发明中的新的非线性滤波算法和 EKF算法计算量的统计结果如表1所示。由于新的非线性滤波算法与EKF算法 的不同仅在于估计误差校正量的求解方程上，为此，表中的结果仅对两种算法 中求解一次估计误差校正量的计算量进行了统计，表中n、m分别表示系统状态 变量的维数和测量输出的维数。

表1 两种算法的计算量统计表

  新的非线性滤波算法   EKF   加法   2n²-2n+mn   2mn²-2mn+2m²n+2n³-n²-n   乘法   2n²+mn   2n³+2mn²+2m²n+mn

EKF滤波增益的计算中存在矩阵求逆，该计算量不便于统计，因而表中对 EKF求解估计误差校正量的计算量统计并未包含矩阵求逆的计算量。由于 m，n≥1，且皆为正整数，对比表中的统计结果可知，EKF的乘法计算量和加法计 算量都分别远大于新的非线性滤波算法的乘法计算量和加法计算量，且随着系 统状态变量维数n和测量输出维数m的增加，EKF算法计算量的增加明显要快 于新的非线性滤波算法的计算量，也即新的非线性滤波算法比EKF在计算的复 杂度上则表现出更明显的优势。由于表1中的对算法计算量的统计结果中未统 计EKF算法中矩阵求逆的计算量，结合前面的两种滤波算法计算量的对比分析 结论可以得到：以算法的计算量与复杂性作为衡量滤波算法的实时性指标时， 新的非线性滤波算法比EKF具有更好的实时性。

在实际仿真中，每一步都是用k时刻的测量去预测k+1时刻的状态估计误 差校正量预测并加入到缓存中，用于计算下一时刻的状态估计误差校正 量，即：此方法能避免缓存上一时刻的测量数据 y_k-1和一步预测

本实施例以状态模型具有典型非线性特性的飞行器再入模型为仿真对象：

${\stackrel{·}{x}}_1\left(t\right)=-x_2\left(t\right)+w_1\left(t\right)$

${\stackrel{·}{x}}_2\left(t\right)=-e^{-\gamma x_1\left(t\right)}{x_2}^2\left(t\right)x_3\left(t\right)+w_2\left(t\right)---\left(10\right)$

${\stackrel{·}{x}}_3\left(t\right)=w_3\left(t\right)$

其中，x₁(t)、x₂(t)、x₃(t)分别表示飞行器的高度、速度及弹道常值系数，w₁(t)、 w₂(t)、w₃(t)为零均值的非相关噪声，其协方差矩阵为Q，γ表示与高度相关的空 气密度常值参数，其值为5×10^-5。雷达的测量模型如下：

$z\left(t\right)=\sqrt{M^2+{[x_1\left(t\right)-H]}^2}+r\left(t\right)---\left(26\right)$

式中，M=10⁵ft表示飞行器与雷达之间的水平距离，H＝10⁵ft表示雷达所在高度， r(t)为零均值的观测噪声，其协方差矩阵为R(t)=10⁴ft²。

系统状态初始值为x(0)=[3×10⁵ 2×10⁴ 10^-3]，为保证仿真中EKF滤波器始 终为无偏估计，仿真中设定系统状态估计的初值与真实的初始值相等，系统离 散化采样时间间隔为0.02s。

由于系统状态模型中的噪声可表征线性化误差，首先考虑飞行器的高度和 速度变化不存在噪声干扰的情况，系统(10)过程噪声的协方差为 Q=diag(0,0,10^-8)。采用本算例中采用新的非线性滤波算法和EKF得到的系统(10) 的三个状态的绝对估计误差对比图如图2~图4所示。图中，横轴表示仿真时间， 图2~图4中的三个子图的纵轴按照从上到下的顺序分别表示新的非线性滤波方 法得到的系统(10)的三个状态的绝对估计误差、EKF方法得到的系统(10)的三个 状态的绝对估计误差和两种方法得到的系统(10)的三个状态的绝对估计误差之 差。由图可知，新的非线性滤波算法得到的三个状态的绝对估计误差随时间收 敛，并逐渐进入稳态，表明本发明中的新的非线性滤波方法是正确的。过渡过 程中两种方法得到的高度和速度的绝对估计误差曲线近乎相同，但进入稳态后， 新的非线性滤波算法的得到的绝对估计误差绝对值比EKF的大。两种方法得到 的弹道常值系数绝对估计误差的数量级皆为10^-8，皆可近似认为无偏估计，但前 20s新的非线性滤波算法得到的弹道常值系数的绝对估计误差比EKF的大，之 后比EKF方法的小，这是因为新的非线性滤波算法中对状态估计校正量的求解 是以高斯噪声到估计误差的传递函数的H2范数最小为设计指标的，对噪声具有 一定的抑制作用，而EKF方法中对x₃无观测，其估计仅依赖于系统模型，过渡 过程中EKF滤波增益为在线调整的，绝对估计误差比新滤波算法的要小，但随 着估计误差逐渐收敛进入稳态后，新的非线性滤波算法对噪声具有抑制作用则 表现出更明显的优势，其绝对估计误差比EKF的小。

以仿真时长130s-200s内的仿真数据为统计对象，以稳态最大绝对估计误差 作为衡量滤波器估计性能的指标，稳态时，新的非线性滤波算法得到的高度、 速度及弹道常系数的最大绝对估计误差值分别为5.393ft，0.0438ft/s、3.979× 10^-8；EK算法得到的高度、速度及弹道常系数的最大绝对估计误差值分别为 2.8374ft，0.0410ft/s、4.845×10^-8，其中新的非线性滤波算法得到的高度与速度 的估计精度近似比EKF分别下降了约0.9倍、0.1倍，但仍满足一般工程上的需 求。

根据表1统计的两种滤波算法的计算量统计结果可知，得到系统(10)三个状 态的一次估计，新的非线性滤波算法和EKF的加法计算量分别为15和60，乘 法计算量分别为21和81，新的非线性滤波算法的计算量远低于EKF。综合本发 明中新的非线性滤波算法和EKF的滤波精度及算法复杂性对比结果可得，本发 明中的新的非线性滤波算法滤波精度虽EKF略差，但由于滤波系数为常值，且 不需要实时计算雅可比矩阵，计算量比EKF大为降低，且算法设计简单，易于 实现，比EKF具有更好的实时性和适用性。