CN103335653A - 火星大气进入段的自适应增量粒子滤波方法 - Google Patents

Abstract

一种火星大气进入段的自适应增量粒子滤波方法，它包括以下步骤：(一)初始化：由先验概率密度函数采样得到粒子并赋予相同的权值；(二)利用各粒子在运动方程和增量量测方程下的似然概率密度函数更新粒子及更新权值；(三)比较步骤二中计算得到的所有权值，得到最大权值、最小权值及相对应的序号并存储；根据增量量测方程和当前时刻对应的增量量测值，求得量测新息计算两最值之间的欧式距离及每个粒子到最小权值所对应粒子的欧式距离；(四)确定自适应系数值并重新计算权值；(五)将步骤四中计算得到的权值进行归一化，得到新的权值；(六)重新采样；(七)返回步骤二，直至时间截止。本方法能消除量测系统中的未知系统误差。

Description

火星大气进入段的自适应增量粒子滤波方法

技术领域

本发明涉及火星大气进入段的自适应增量粒子滤波方法。属于航天导航技术领域。

背景技术

粒子滤波方法是非常常见的一种确定航天器位置速度方法。其中测量方法中未知测量系统误差的大小对航天器位置速度最终状态的确定影响很大，过大的未知测量系统误差数据会导致位置速度误差的增大，增大航天器的导航误差、降低导航精度。

现有技术中，可以用于确定航天器位置速度的方法有多种。

现有技术一，基于泰勒展开的扩展卡尔曼滤波估计方法。该方法给出了非线性动力学方程和非线性测量方程的泰勒展开加权融合的估算公式。

现有技术二，基于sigma点集(为正态分布采样策略)的无迹卡尔曼滤波方法。先根据正态分布的均值和方差计算出sigma点集，并确定出各点的权值，再通过动力学方程计算出航天器的位置速度，然后通过量测方程得到的量测数据对航天器的位置速度进行调整修正。

现有技术一在超音速强耦合强干扰非线性环境中将动力学展开得到显著的误差，因此不太适用于火星大气进入段。

现有技术二在测量手段有限，测量数据少，难以实现系统误差的估算和消除。测量设备即使在地面不同环境中产生的系统误差各不相同。即使在地面试验中已经校准的系统误差在新的火星环境(其环境不同与地球)中不再准确，因此量测数据中的系统误差将影响航天器的位置速度进行调整修正，因此不太适用于火星大气进入段。

现有技术一二对于火星环境中未知的，有偏的和非正态的噪声统计量得到显著的误差，因此不太适用于火星大气进入段。

发明内容

1、目的：本发明的目的是提供一种火星大气进入段的自适应增量粒子滤波方法，以减小航天器位置速度误差，提高其精度。

2、技术方案：本发明的目的是通过以下技术方案来实现的。

工程实践中所遇到状态方程和量测方程往往是非线性的，为了描述系统状态估计问题，定义非线性离散增量系统如下：

运动方程        x_k=f_k-1(x_k-1,w_k-1)

传统量测方程    z_k=h_k(x_k)+b_k+v_k

增量量测方程    Δz_k=h_k(x_k,x_k-1)+v_k

式中,x_k为状态向量；f_k(·)和h_k(·)为非线性向量离散函数；w_k为系统噪声向量；Δz_k=z_k-z_k-1为量测向量增量；z_k为量测向量；v_k为量测噪声向量；b_k为未知的量测系统误差。两次相邻量测时，量测的系统误差相同或相近，两量测值之差对应的量测系统误差为零或较小量。本滤波方法基于增量量测方程进行自适应增量粒子滤波。本滤波方法既能适应环境变化又能消除量测过程中的系统误差，最终提高滤波精度。

本发明一种火星大气进入段的自适应增量粒子滤波方法，它包括以下步骤：

步骤一、初始化：时刻k=0,由先验概率密度函数p(x₀)采样,得到粒子并赋予相同的权值 ${\mathrm{\ω}}_0^i=\frac{1}{N},i=1,...,N;$

式中：k为采样步，N为采样粒子总数，p(x₀)为采样的先验概率密度函数，如正态分布，伽马分布，卡方分布等等；

步骤二、更新粒子及更新权值：

(1)从

中随机抽取N个有限粒子；

(2)逐点计算对应的 $p\left(x_k^i\right|x_{k-1}^i)$ 和 $p\left({\mathrm{\Δz}}_k\right|x_k^i)$ ，即 $x_k^{i^{*}}=p\left(x_k\right|x_{k-1}^i)i=1,...,N,$ 计算每个粒子的似然概率密度函数

(3)利用下式计算对应粒子的重要性权值：

${\mathrm{\ω}}_k^i={\mathrm{\ω}}_{k-1}^i\frac{p\left({\mathrm{\Δz}}_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,{\mathrm{\Δz}}_k)},=1,...,N,$ 即更新权值 ${\mathrm{\ω}}_k^{i^{*}}=p\left({\mathrm{\Δz}}_k\right|x_k^{i^{*}});$

式中：利用带有权值系数的粒子

来描述k时刻的验后概率密度p(x_0:k|Δz_1:k)；由于p(x_0:k|Δz_1:k)采样困难,故可以使用重要性函数q(x_0:k|Δz_1:k)进行采样；粒子可由重要性函数q(x_0:k|Δz_1:k)得到，则权值为

若重要性函数q(x_0:k|Δz_1:k)可以分解为q(x_0:k|Δz_1:k)=q(x_k|x_0:k-1,Δz_1:k)q(x_0:k-1|Δz_1:k-1)，并且p(x_0:k|Δz_1:k)可分解为p(x_0:k|Δz_1:k)=p(Δz_1:k|x_k)p(x_k|x_k-1)×p(x_0:k-1|Δz_1:k-1)则重要性权值更新公式为

${\mathrm{\ω}}_k^i\∝{\mathrm{\ω}}_{k-1}^i\frac{p\left({\mathrm{\Δz}}_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,{\mathrm{\Δz}}_k)};$

为k时刻各粒子的概率形式，即运动方程在离散时间对应的Markov随机过程，由动力学方程f(·)和过程噪声分布p(v_k)决定；

为k时刻各粒子在增量量测方程下的似然概率密度函数，由量测系统h(·)和量测噪声分布p(w_k)决定。

q(x_k|x_0:k-1,Δz_1:k)为重要性分布函数；

p(x_0:k|Δz_1:k)为状态更新方程的验后概率密度。

步骤三、比较步骤二中计算得到的所有权值，可以得到最大权值、最小权值

和相对应的序号i_ω-max,i_ω-min并存储；根据增量量测方程Δz_k=h_k(x_k,x_k-1)+v_k和当前时刻对应的增量量测值Δz_k，求得量测新息 ${\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\max }}={\mathrm{\Δz}}_k-h\left(x_k^{i_{\mathrm{\ω}-\max }^{*}}{\hat{x}}_{k-1}\right),{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\min }}={\mathrm{\Δz}}_k-h\left(x_k^{i_{\mathrm{\ω}-\min }^{*}}{\hat{x}}_{k-1}\right)$ 计算与

之间的欧式距离L_max及每个粒子到最小权值所对应粒子的欧式距离Li；

式中：

为k-1时刻的状态估计值；

欧式距离分别按照以下两个式子进行计算：

$L_{\max }=({\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\max }}-{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\min }}){({\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\max }}-{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\min }})}^T,$

$L^i=({\mathrm{\Δz}}_k-{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}\min }}){({\mathrm{\Δz}}_k-{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}\min }})}^T;$

步骤四、计算确定自适应系数β值，由量测噪声的统计特性决定；在自适应算法中β值尤为重要；当量测噪声较低时，β=0，不对似然分布作调整；当量测噪声较高，即当似然分布呈尖峰状态或位于转移先验分布尾部时，β＞0，即人为使得似然分布较广一些； $\mathrm{\β}=\left\{\begin{array}{cc}K/\mathrm{\α}& \mathrm{\α}\≤\mathrm{\ϵ}\\ 0& \mathrm{\α}>\mathrm{\ϵ}\end{array}\right.$ 式中，ε为阀值，根据经验确定；K为比例常数，K/α＞0。在本算法中取β=1；然后重新计算权值，其方式如下：

${\mathrm{\ω}}_k^{i^{*}}={\mathrm{\ω}}_k^{i^{*}}+({\mathrm{\ω}}_k^{i_{\mathrm{\ω}-\max }^{*}}/N)\·\sin (L^i/L_{\max }\·\mathrm{\π}/2)\·\mathrm{\β};$

步骤五、将步骤四中计算得到的权值进行归一化，得到新的权值，其方法如下：

${\mathrm{\ω}}_k^{i^{*}}={\mathrm{\ω}}_k^{i^{*}}/\underset{j}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{j^{*}},$

将步骤四中计算得到的每一个权值除以步骤四中计算得到的所有权值的总和得到归一化后新的权值，然后利用步骤五计算得到的权值与对应的粒子相乘再求和，则可得k时刻未知参数x_k的最小均方估计

步骤六、重新采样：对

进行重采样得到新粒子群

及其权值

步骤七、令k=k+1，返回步骤二往下进行，直到k等于火星大气进入时间截止对应的时刻T时；至此完成火星大气进入段的自适应增量粒子滤波方法。

其中，在步骤一中所述的先验概率一般情况下概率密度函数为正态分布函数，具体的均值和方差由火星实际工程背景估算得到；其中的粒子(即随机样本)数理论上是越大越能反映实际情况，但是火星着陆时间短，只能选取恰当的数据量既能反映真实的情况，也能使得计算时间满足火星着陆要求。

其中，在步骤二中采用序贯重要性采样(Sequential importance sampling,SIS)方法来构造未知系统状态的概率密度函数；即利用系统生成的粒子(即随机样本)的加权组合来表示相应的验后概率密度函数，并利用这些粒子(即随机样本)和权值得到状态的估计值；已知系统的状态先验条件概率p(x₀),利用带有权值系数的粒子来描述k时刻的验后概率密度p(x_0:k|Δz_1:k),则k时刻的验后概率密度可离散地加权为 $p\left(x_{0:k}\right|{\mathrm{\Δz}}_{1:k})\≈\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i\mathrm{\δ}(x_{0:k}-x_{0:k}^i),$ 式中,权值系数 ${\mathrm{\ω}}_k^i$ 是相应的随机粒子 $x_k^i$ 的归一化权值,即由于p(x_0:k|Δz_1:k)采样困难,故可以使用重要性函数q(x_0:k|Δz_1:k)进行采样；若粒子

可由q(x_0:k|Δz_1:k)得到,则权值为

若q(x_0:k|Δz_1:k)可以分解为q(x_0:k|Δz_1:k)=q(x_k|x_0:k-1,Δz_1:k)q(x_0:k-1|Δz_1:k-1)，并且p(x_0:k|Δz_1:k)可分解为p(x_0:k|Δz_1:k)=p(Δz_1:k|x_k)p(x_k|x_k-1)×p(x_0:k-1|Δz_1:k-1)则重要性权值值更新公式为 ${\mathrm{\ω}}_k^i\∝{\mathrm{\ω}}_{k-1}^i\frac{p\left({\mathrm{\Δz}}_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,{\mathrm{\Δz}}_k)}$ 从而,验后概率密度的加权近似为 $p\left(x_k\right|{\mathrm{\Δz}}_{1:k})\≈\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i\mathrm{\δ}(x_k-x_k^i).$

其中，在步骤三中计算的欧式距离为最大权值粒子和最小权值粒子之间的距离，他们是采样后的粒子的上界和下界，即确定了粒子大小的范围及概率特征，通过改变粒子对应权值的分布得到改变似然函数的分布自适应增量粒子滤波。

其中，在步骤四中的自适应系数β，由量测噪声的统计特性决定；在自适应算法中β值尤为重要；当量测噪声较低时，β=0，不对似然分布作调整；当量测噪声较高，即当似然分布呈尖峰状态或位于转移先验分布尾部时，β＞0，即人为使得似然分布较广一些； $\mathrm{\β}=\left\{\begin{array}{cc}K/\mathrm{\α}& \mathrm{\α}\≤\mathrm{\ϵ}\\ 0& \mathrm{\α}>\mathrm{\ϵ}\end{array}\right.$ 式中，ε为阀值，根据经验确定；K为比例常数，K/α＞0。在本算法中取β=1。

3、优点和功效：

本发明统筹考虑了火星实际大气进入过程中，非线性、非高斯随机系统在量测系统存在未知系统误差条件下的航天器位置速度估计问题。通过自适应改变粒子对应权值的分布得到改变似然函数的分布自适应增量粒子滤波，避免粒子退化实现正确表达实际的验后概率分布。并在计算过程中引入了测量增量值，消除了量测系统中的未知系统误差。因而本发明提出的算法可以有效保证航天器在火星大气进入段的位置速度估计。

附图说明

图1为非线性增量系统的概率密度函数图解模型

图2为真实状态、跟踪估计值

图3为跟踪估计值与真实状态值之间的误差比较

图4为本发明自适应增量粒子滤波方法流程图

图中的代号、符号说明如下：

Δz为量测增量向量，Δz_k=z_k-z_k-1，z_k为量测向量。

x_k为状态向量。

p(x_k|x_k-1)为该火星大气进入非线性系统的马尔科夫过程概率密度分布函数，由动力学方程f(·)和系统噪声分布p(v_k)决定。

p(Δz_k|x_k)为火星大气进入量测向量增量的概率密度分布函数，由量测方程h(·)和量测噪声分布p(w_k)决定。

True state为状态真实值。

AIPF为自适应增量粒子滤波方法。

APF为自适应粒子滤波方法。

具体实施方式

本发明涉及火星大气进入段的自适应增量粒子滤波方法，具体实施步骤如下：

航天器沿飞行轨迹进入火星大气在飞行力学按照火星环境设定其动力学系统参数和系统噪声分布。而量测系统是根据火星大气进入实际工程中的仅有的测量方式建立。其量测噪声也有测量设备分析估计近似得到，系统误差在测量数据较少的情况下很难确定成为未知量测系统误差。动力学系统对应的非线性运动学方程和量测方程如下：

运动方程 $x_k=10+0.5x_{k-1}+\frac{{2.5x}_{k-1}}{{1+x}_{k-1}^2}+2\cos \left(1.2(k-1)\right)+w_{k-1}$

量测方程 $z_k=\frac{x_k^2}{5}+a{+v}_k$

增量量测方程    Δz_k=z_k-z_k-1+v_k

式中，w_k和v_k是独立的高斯白噪声序列，且量测噪声v_k的均值r=0，R=1；量测系统误差a=3为未知量，采样粒子数n=150.系统噪声w_k的均值q=0和方差Q的变化如下式

$Q=\left\{\begin{array}{cc}0.1& k\∈\left[\mathrm{1,50}\right)\\ 0.5& k\∈\left[\mathrm{50,100}\right)\\ 0.2& k\∈\left[\mathrm{100,150}\right]\end{array}\right.$

在时间截止对应的时刻T之前，本滤波方法按照图1非线性增量系统的概率密度函数图解模型顺序进行计算，每个时刻按照图4自适应增量粒子滤波具体实施的流程图进行滤波。

步骤一：设定粒子数N=200，根据各个状态量的噪声分布设定的先验概率p(x₀)～N(0,2)采样，得到200个粒子:i=1,...,200，并赋予相同的权值

粒子数N越大，反映实际状态越准确，但是计算量较大。在计算速度和有限时间容许的范围内，可以进一步提高粒子数量。

步骤二、更新粒子及更新权值：

中随机抽取N个有限粒子；根据运动方程f(·)进行状态更新；

(2)逐点计算N个粒子对应的增量量测估计

(3)利用下式计算对应粒子的重要性权值：

${\mathrm{\ω}}_k^i=\left(\frac{\exp (-\frac{{\mathrm{\Δz}}_k-{\widehat{\mathrm{\Δz}}}_k^i}{2^{*}2})}{\sqrt{2^{*}2\mathrm{\π}}}\right),i=1,...,N,$ 即更新权值 ${\mathrm{\ω}}_k^{i^{*}}=p\left({\mathrm{\Δz}}_k\right|x_k^{i^{*}});$

步骤三、比较步骤二中计算得到的所有权值，可以得到最大权值

、最小权值

和相对应的序号i_ω-max,i_ω-min并存储。根据增量量测方程Δz_k=h_k(x_k,x_k-1)+v_k和当前时刻对应的增量量测值Δz_k，求得量测新息 ${\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\max }}={\mathrm{\Δz}}_k-h\left(x_k^{i_{\mathrm{\ω}-\max }^{*}}{\hat{x}}_{k-1}\right),{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\min }}={\mathrm{\Δz}}_k-h\left(x_k^{i_{\mathrm{\ω}-\min }^{*}}{\hat{x}}_{k-1}\right),$ 计算

与

之间的欧式距离L_max及每个粒子到最小权值所对应粒子的欧式距离Lⁱ。

欧式距离分别按照以下两个式子进行计算：

$L_{\max }=({\mathrm{\Δ}}_k^{i_{\mathrm{\ω}\max }}-{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\min }}){\left({\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\max }}{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}\min }}\right)}^T,$

$L^i=({\mathrm{\Δz}}_k-{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}\min }}){({\mathrm{\Δz}}_k-{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}\min }})}^T.$

步骤四：由于量测噪声较高，自适应系数β值取1进行计算。然后按照公式 ${\mathrm{\ω}}_k^{i^{*}}={\mathrm{\ω}}_k^{i^{*}}+({\mathrm{\ω}}_k^{i_{\mathrm{\ω}-\max }^{*}}/N)\·\sin (L^i/L_{\max }\·\mathrm{\π}/2)$ 重新进行计算权值。

步骤五：，利用

将步骤四中重新计算得到的权值进行归一化处理，得到新的权值并与对应的粒子计算出状态估计量。

步骤六、重采样。随机生成一个0～1之间的随机数，将步骤五中得到的新权值从第一个开始与其比较：如果大于等于随机数就取相对应的粒子，如果小于随机数就将第一个权值往后累加直到累加得到的权值和大于等于随机数就取此时相对应的粒子。共进行N次。重采样结束。

步骤七、令k=k+1，返回步骤二往下进行。直到k等于截止对应的时刻T时。至此完成自适应增量粒子滤波方法。

从图2中可以看到，自适应增量粒子滤波方法得到的状态估计比自适应粒子滤波方法的精度高，能够消除系统误差对滤波的影响，提高滤波的稳定性。图3给出两种滤波结果与真值的误差比较，可以明显看到，自适应增量粒子滤波方法得到的状态估计比自适应粒子滤波方法的误差要小得多，能很好跟踪状态x_k的变化。

以上所述仅为本发明较佳的实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化和替换都应涵盖在本发明的保护范围之内，另外本发明提供的方法可以集成到火星大气进入航天器位置速度估计软件中。

Claims (5)

1.一种火星大气进入段的自适应增量粒子滤波方法，其特征在于：它包括以下步骤：

步骤一、初始化：时刻k=0,由先验概率密度函数p(x₀)采样,得到粒子

:i=1,…,N，并赋予相同的权值

式中：k为采样步，N为采样粒子总数，p(x₀)为采样的先验概率密度函数，如正态分布，伽马分布，卡方分布；

步骤二、更新粒子及更新权值：

(1)从

中随机抽取N个有限粒子；

(2)逐点计算对应的 $p\left(x_k^i\right|x_{k-1}^i)$ 和 $p\left({\mathrm{\Δz}}_k\right|x_k^i),$ 即 $x_k^{i^{*}}=p\left(x_k\right|x_{k-1}^i),i=1,...,N,$ 计算每个粒子的似然概率密度函数

(3)利用下式计算对应粒子的重要性权值：

${\mathrm{\ω}}_k^i={\mathrm{\ω}}_{k-1}^i\frac{p\left({\mathrm{\Δz}}_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,{\mathrm{\Δz}}_k)},i=1,...,N,$ 即更新权值 ${\mathrm{\ω}}_k^{i^{*}}=p\left({\mathrm{\Δz}}_k\right|x_k^{i^{*}});$

式中：利用带有权值系数的粒子

来描述k时刻的验后概率密度p(x_0:k|Δz_1:k)；由于p(x_0:k|Δz_1:k)采样困难,故可以使用重要性函数q(x_0:k|Δz_1:k)进行采样；粒子

可由重要性函数q(x_0:k|Δz_1:k)得到，则权值为；若重要性函数q(x_0:k|Δz_1:k)分解为q(x_0:k|Δz_1:k)=q(x_k|x_0:k-1,Δz_1:k)q(x_0:k-1|Δz_1:k-1)，并且p(x_0:k|Δz_1:k)分解为p(x_0:k|Δz_1:k)=p(Δz_1:k|x_k)p(x_k|x_k-1)×p(x_0:k-1|Δz_1:k-1)则重要性权值更新公式为

${\mathrm{\ω}}_k^i\∝{\mathrm{\ω}}_{k-1}^i\frac{p\left({\mathrm{\Δz}}_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,{\mathrm{\Δz}}_k)};$

为k时刻各粒子的概率形式，即运动方程在离散时间对应的Markov随机过程，由动力学方程f(·)和过程噪声分布p(v_k)决定；

为k时刻各粒子在增量量测方程下的似然概率密度函数，由量测系统h(·)和量测噪声分布p(w_k)决定；

q(x_k|x_0:k-1,Δz_1:k)为重要性分布函数；

p(x_0:k|Δz_1:k)为状态更新方程的验后概率密度；

步骤三、比较步骤二中计算得到的所有权值，得到最大权值

、最小权值

和相对应的序号i_ω-max,i_ω-min并存储；根据增量量测方程Δz_k=h_k(x_k,x_k-1)+v_k和当前时刻对应的增量量测值Δz_k，求得量测新息

${\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\min }}={\mathrm{\Δz}}_k-h(x_k^{i_{\mathrm{\ω}-\min }^{*}},{\hat{x}}_{k-1}),$ 计算

与

之间的欧式距离L_max及每个粒子到最小权值所对应粒子的欧式距离Li；

式中：为k-1时刻的状态估计值；

欧式距离分别按照以下两个式子进行计算：

$L_{\max }=({\mathrm{\Δ}}_k^{i_{\mathrm{\ω}\max }}-{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\min }}){\left({\mathrm{\Δz}}_k^{i_{\mathrm{\ω}-\max }}{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}\min }}\right)}^T,$

$L^i=({\mathrm{\Δz}}_k-{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}\min }}){({\mathrm{\Δz}}_k-{\mathrm{\Δz}}_k^{i_{\mathrm{\ω}\min }})}^T;$

步骤四、计算确定自适应系数β值，由量测噪声的统计特性决定；在自适应算法中β值尤为重要；当量测噪声较低时，β=0，不对似然分布作调整；当量测噪声较高，即当似然分布呈尖峰状态或位于转移先验分布尾部时，β＞0，即人为使得似然分布较广一些； $\mathrm{\β}=\left\{\begin{array}{cc}K/\mathrm{\α}& \mathrm{\α}\≤\mathrm{\ϵ}\\ 0& \mathrm{\α}>\mathrm{\ϵ}\end{array}\right.$ 式中，ε为阀值，根据经验确定；K为比例常数，K/α＞0；在本算法中取β=1；然后重新计算权值，其方式如下：

${\mathrm{\ω}}_k^{i^{*}}={\mathrm{\ω}}_k^{i^{*}}+({\mathrm{\ω}}_k^{i_{\mathrm{\ω}-\max }^{*}}/N)\·\sin (L^i/L_{\max }\·\mathrm{\π}/2)\·\mathrm{\β};$

步骤五、将步骤四中计算得到的权值进行归一化，得到新的权值，其方法如下：

${\mathrm{\ω}}_k^{i^{*}}={\mathrm{\ω}}_k^{i^{*}}/\underset{j}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{j^{*}},$

将步骤四中计算得到的每一个权值除以步骤四中计算得到的所有权值的总和得到归一化后新的权值，然后利用步骤五计算得到的权值与对应的粒子相乘再求和，则得到k时刻未知参数x_k的最小均方估计

步骤六、重新采样：对

进行重采样得到新粒子群 $\{x_k^i,{\mathrm{\ω}}_k^i:i=1,...,N\},$ 及其权值 ${\mathrm{\ω}}_k^i=1/N;$

步骤七、令k=k+1，返回步骤二往下进行，直到k等于火星大气进入时间截止对应的时刻T时；至此完成火星大气进入段的自适应增量粒子滤波方法。

2.根据权利要求1所示的一种火星大气进入段的自适应增量粒子滤波方法，其特征在于：在步骤一中所述的先验概率一般情况下概率密度函数为正态分布函数，具体的均值和方差由火星实际工程背景估算得到；其中的粒子即随机样本数理论上是越大越能反映实际情况，但是火星着陆时间短，只能选取恰当的数据量既能反映真实的情况，也能使得计算时间满足火星着陆要求。

3.根据权利要求1所示的一种火星大气进入段的自适应增量粒子滤波方法，其特征在于：在步骤二中采用序贯重要性采样即SIS方法来构造未知系统状态的概率密度函数；利用系统生成的粒子即随机样本的加权组合来表示相应的验后概率密度函数，并利用这些粒子即随机样本和权值得到状态的估计值；已知系统的状态先验条件概率p(x₀),利用带有权值系数的粒子来描述k时刻的验后概率密度p(x_0:k|Δz_1:k),则k时刻的验后概率密度离散地加权为 $p\left(x_{0:k}\right|{\mathrm{\Δz}}_{1:k})\≈\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i\mathrm{\δ}(x_{0:k}-x_{0:k}^i),$ 式中,权值系数是相应的随机粒子

的归一化权值,即

由于p(x_0:k|Δz_1:k)采样困难,故可以使用重要性函数q(x_0:k|Δz_1:k)进行采样；若粒子

由q(x_0:k|Δz_1:k)得到,则权值为

若q(x_0:k|Δz_1:k)分解为q(x_0:k|Δz_1:k)=q(x_k|x_0:k-1,Δz_1:k)q(x_0:k-1|Δz_1:k-1)，并且p(x_0:k|Δz_1:k)分解为p(x_0:k|Δz_1:k)=p(Δz_1:k|x_k)p(x_k|x_k-1)×p(x_0:k-1|Δz_1:k-1)则重要性权值更新公式为 ${\mathrm{\ω}}_k^i\∝{\mathrm{\ω}}_{k-1}^i\frac{p\left({\mathrm{\Δz}}_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,{\mathrm{\Δz}}_k)}$ 从而,验后概率密度的加权近似为

$p\left(x_k\right|{\mathrm{\Δz}}_{1:k})\≈\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^i\mathrm{\δ}(x_k-x_k^i).$

4.根据权利要求1所示的一种火星大气进入段的自适应增量粒子滤波方法，其特征在于：在步骤三中计算的欧式距离为最大权值粒子和最小权值粒子之间的距离，他们是采样后的粒子的上界和下界，即确定了粒子大小的范围及概率特征，通过改变粒子对应权值的分布得到改变似然函数的分布自适应增量粒子滤波。

5.根据权利要求1所示的一种火星大气进入段的自适应增量粒子滤波方法，其特征在于：在步骤四中所述的自适应系数β，由量测噪声的统计特性决定；在自适应算法中β值尤为重要；当量测噪声较低时，β=0，不对似然分布作调整；当量测噪声较高，即当似然分布呈尖峰状态或位于转移先验分布尾部时，β>0，即人为使得似然分布较广一些； $\mathrm{\β}=\left\{\begin{array}{cc}K/\mathrm{\α}& \mathrm{\α}\≤\mathrm{\ϵ}\\ 0& \mathrm{\α}>\mathrm{\ϵ}\end{array}\right.$ 式中，ε为阀值，根据经验确定；K为比例常数，K/α＞0；在本算法中取β=1。

Images