行星大气进入着陆器落点不确定度分析方法

摘要

本发明属于航天器着陆与返回技术领域，涉及一种行星大气进入着陆器落点不确定度分析方法。首先根据系统初始状态的不确定性分布将状态用Askey正交多项式逼近，然后将状态带入到系统动力学中，根据Galerkin投影法则，将表示原系统的随机微分方程转化为一个等效的高维确定性微分方程，最后利用龙格-库塔等数值积分方法，得到各时刻表示系统状态的正交多项式系数，从而得到系统状态的统计特性，并且在整个过程中根据着陆器状态的统计特性自适应调整正交多项式基底，克服截断误差带来的影响。该发明能够准确的估计系统状态的统计特性，并且计算效率明显提高。

说明书

技术领域

本发明属于航天器着陆与返回技术领域，涉及一种行星大气进入着陆器落点不确定度分析方法。

背景技术

在带有大气行星上完成着陆任务，需要在任务前选定预定的着陆点，但探测器在火星大气进入点处的导航控制误差，探测器的气动参数以及火星大气模型的不确定性，都会严重影响着陆器最终的着陆精度，甚至关乎任务的成败。因此，分析这些偏差以及不确定性对着陆点的影响，是一项必不可少的工作；针对带有大气的行星着陆任务，发展一种快速的落点不确定度分析方法，对降低未来火星着陆设计周期和成本，提高设计效率很有意义。目前，在处理这个问题的方法中，总的来说有三类，一是根据系统状态初值及系统方程中不确定参数的统计特性，选择足够多的采样点，进行蒙特卡洛仿真，从而得到各个时刻系统状态的统计特性；二是将系统方程进行线性化，利用线性系统理论对着陆点的统计特性进行分析；三是利用根据系统初始状态的不确定性分布将状态用Askey正交多项式逼近，然后将状态带入到系统动力学中，根据Galerkin投影法则，将表示原系统的随机微分方程转化为一个等效的高维确定性微分方程，最后利用龙格-库塔等数值积分方法，得到各时刻表示系统状态的正交多项式系数，从而得到系统状态的统计特性。

第一类方法需要较高的计算代价，利用这类方法往往需要较长的任务周期，第二类方法虽然计算效率高，但线性化地方法使得在系统初始状态偏差较大时，出现发散现象；第三类方法有完整的数学理论体系，并且计算效率比较高，具有进一步发展的潜力。参见AvinashPrabhakar，James Fisher and Raktim Bhattacharya.Polynomial Chaos-Based Analysis ofProbabilistic Ucertainty in Hypersonic Flight Dynamics[J].Journal of Guidance，Control，and Dynamics.2010，33(1)：222-234.中，利用Askey正交多项式和Galerkin投影法将系统表示为等价的高阶微分方程来求解着陆器状态的统计特性，但其没有考虑用Askey正交多项式表示着陆器状态时的截断误差，从而导致在多误差源影响下，算法容易发散的问题。

发明内容

本发明针对现有的行星大气进入着陆器落点不确定性分析技术存在的计算效率低的情况，提出一种行星大气进入着陆器落点不确定度分析方法，能够准确的估计系统状态的统计特性，并且计算效率明显提高。

该行星大气进入着陆器落点不确定度分析方法：

第一步：根据系统初始状态的不确定性分布将状态用Askey正交多项式逼近，构建正交多项式基；

第二步：将系统状态和不确定参数带入到系统动力学中，将表示原系统的随机微分方程转化为一个等价的高阶确定性微分方程；

第三步：利用龙格-库塔等数值积分算法对此高阶确定性微分方程进行积分，求解确定性微分方程，同时对求得的逼近着陆器状态的正交多项式系数进行检测，如果正交多项式的非线性项系数超过预定比例，那么进入第四步，否则进入第五步；

第四步：根据此时的着陆器状态分布特性，利用施密特正交化办法构建新的正交多项式，用新的正交多项式逼近此时的着陆器状态，从新转化成等效确定性微分方程，利用龙格-库塔方法对其进行积分，并监测非线性项系数与线性项系数的比例；

第五步：利用施密特正交化办法建立新的正交多项式，以此类推，直至所需要的停止条件；

第六步：利用数学期望和数学方差的定义，结合每个时刻表示状态的正交多项式，计算此时系统状态的统计特性。

本发明的有益效果：

该发明针对大气进入类行星着陆器落点不确定度问题，能够确保对着陆器统计特性的快速准确估计，并且克服了在多误差源干扰情况下算法发散的问题。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合附图对本发明的实施例作详细说明：本实施例在以本发明的技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

本部分以火星着陆落点偏差的不确定度分析问题为例，给出具体的实施方式。

火星着陆系统动力学为：

(1)

$>+v\cos \gamma /(R_m+h)$>

其中，h表示着陆器距离火星表面的距离，v表示着陆器速度的大小，γ表示航迹角，μ表示火星引力系数，R_m表示火星半径，B表示着陆器的弹道系数，k表示着陆器的升阻比，φ表示倾侧角，λ表示大气模型不确定性因子，ρ表示火星大气密度，其与着陆器距离火星表面高度的关系如式(2)所示，它是根据NASA开发的火星大气模型MarsGram所生成的数据进行最小二乘拟合得到的。

T＝1.4×10^-13h³-8.85×10^-9h²

-1.245×10^-3h+205.3645                     (2)

P＝559.351005946503e^-0.000105h

ρ＝P/188.95110711075T

假设系统初始状态及不确定性参数的标称状态及不确定性如下表所示

则本系统状态在300s内统计特性可以按照以下方式求取：

步骤1：根据高斯分布的概率密度函数和施密特正交化算法，构建正交多项式基Hi；

步骤2：将系统状态和不确定参数表示成以下形式，

B(ζ)＝B₀H₀(ζ)+B₁H₁(ζ)；k(ζ)＝k₀H₀(ζ)+k₁H₁(ζ)

$>\lambda \left(\zeta \right)={\lambda }_0{H}_0\left(\zeta \right)+{\lambda }_1{H}_1\left(\zeta \right);h(t,\zeta )=\underset{0}{\overset{P}{\Sigma }}{h}_t\left(t\right)H_t\left(\zeta \right)$>

$>v(t,\zeta )=\underset{0}{\overset{P}{\Sigma }}{v}_i\left(t\right)H_i\left(\zeta \right);$>$>\gamma (t,\zeta )=\underset{0}{\overset{P}{\Sigma }}{\gamma }_i\left(t\right)H_i\left(\zeta \right)$>

步骤3：利用龙格-库塔方法对下式进行积分

$>\frac{1}{⟨H_m^2⟩}⟨\cos \left(\underset{i=0}{\overset{P}{\Sigma }}{\gamma }_i{H}_i\right)\underset{i,j=0}{\overset{P}{\Sigma }}{z}_i{v}_j{H}_i{H}_j{H}_m⟩$>

步骤四：若在积分过程中满足

max(|h₂(t₁)|，L，|h_p(t₁)|)≥|h₁(t₁)|/θ₁    or

max(|v₂(t₁)|，L，|v_p(t₁)|)≥|v₁(t₁)|/θ₂    or

max(|γ₂(t₁)|，L，|γ_p(t₁)|)≥|γ₁(t₁)|/θ₃

则转第五步，若不满足，则返回步骤三；

步骤五：根据施密特正交化办法，及下式所表示状态的统计特性，构建新的正交基ξ₁，ξ₂，ξ₃；

$>{\xi }_1=\underset{0}{\overset{P}{\Sigma }}{h}_i\left(t_1\right)H_i\left(\zeta \right)=T_1\left(\zeta \right)$>

$>{\xi }_2=\underset{0}{\overset{P}{\Sigma }}{v}_i\left(t_1\right)H_i\left(\zeta \right)=T_2\left(\zeta \right)$>

$>{\xi }_3=\underset{0}{\overset{P}{\Sigma }}{\gamma }_i\left(t_1\right)H_i\left(\zeta \right)=T_3\left(\zeta \right)$>

步骤六：构建新的正交多项式，

$>h(t,{\xi }_1,{\xi }_2,{\xi }_3)=\underset{0\le l+m+n\le P}{\overset{P}{\Sigma }}{h}_{\mathrm{lmn}}\left(t\right){\phi }_l^{{\xi }_1}\left({\xi }_1\right){\phi }_m^{{\xi }_2}\left({\xi }_2\right){\phi }_n^{{\xi }_3}\left({\xi }_3\right)$>

$>v(t,{\xi }_1,{\xi }_2,{\xi }_3)=\underset{0\le l+m+n\le P}{\overset{P}{\Sigma }}{v}_{\mathrm{lmn}}\left(t\right){\phi }_l^{{\xi }_1}\left({\xi }_1\right){\phi }_m^{{\xi }_2}\left({\xi }_2\right){\phi }_n^{{\xi }_3}\left({\xi }_3\right)$>

$>\gamma (t,{\xi }_1,{\xi }_2,{\xi }_3)=\underset{0\le l+m+n\le P}{\overset{P}{\Sigma }}{\gamma }_{\mathrm{lmn}}\left(t\right){\phi }_l^{{\xi }_1}\left({\xi }_1\right){\phi }_m^{{\xi }_2}\left({\xi }_2\right){\phi }_n^{{\xi }_3}\left({\xi }_3\right)$>

步骤七：对新建立的多项式赋初值

$>h_{\mathrm{lmn}}\left(t_1\right)=\left(\begin{array}{cc}{-\phi }_0^{\left({\xi }_1\right)}& \mathrm{ifl}=m=n=0,\\ 1& \mathrm{ifl}=1^m=n=0,\\ 0& \mathrm{otherwise},\end{array}\right)$>

$>v_{\mathrm{lmn}}\left(t_1\right)=\left(\begin{array}{cc}{-\phi }_0^{\left({\xi }_2\right)}& \mathrm{ifl}=m=n=0,\\ 1& \mathrm{ifm}=1^l=n=0,\\ 0& \mathrm{otherwise},\end{array}\right)$>

$>{\gamma }_{\mathrm{lmn}}\left(t_1\right)=\left(\begin{array}{cc}{-\phi }_0^{\left({\xi }_3\right)}& \mathrm{ifl}=m=n=0,\\ 1& \mathrm{ifn}=1^l=n=0,\\ 0& \mathrm{otherwise},\end{array}\right)$>

步骤八：对式

$>\frac{\mu }{{(R_m+\underset{0\le l+m+n\le P}{\Sigma }{v}_{\mathrm{lmn}}{\phi }_l^{\left({\xi }_1\right)}{\phi }_m^{\left({\xi }_2\right)}{\phi }_n^{\left({\xi }_3\right)})}^2\times \sin \left(\underset{0\le l+m+n\le P}{\Sigma }{\gamma }_{\mathrm{lmn}}{\phi }_l^{\left({\xi }_1\right)}{\phi }_m^{\left({\xi }_2\right)}{\phi }_n^{\left({\xi }_{\text{3}}\right)}\right)}$>

$>\frac{\mu \cos \left(\underset{0\le l+m+n\le P}{\Sigma }{\gamma }_{\mathrm{lmn}}{\phi }_l^{\left({\xi }_1\right)}{\phi }_m^{\left({\xi }_2\right)}{\phi }_n^{\left({\xi }_3\right)}\right)}{\underset{0\le l+m+n\le P}{\Sigma }{v}_{\mathrm{lmn}}{\phi }_l^{\left({\xi }_1\right)}{\phi }_m^{\left({\xi }_2\right)}{\phi }_n^{\left({\xi }_3\right)}{(R_m+\underset{0\le l+m+n\le P}{\Sigma }{v}_{\mathrm{lmn}}{\phi }_l^{\left({\xi }_1\right)}{\phi }_m^{\left({\xi }_2\right)}{\phi }_n^{\left({\xi }_3\right)})}^2}+$>

$>\frac{\cos \mathrm{\phi k\lambda \rho }}{2B}\times \underset{0\le l+m+n\le P}{\Sigma }{v}_{\mathrm{lmn}}{\phi }_l^{\left({\xi }_1\right)}{\phi }_m^{\left({\xi }_2\right)}{\phi }_n^{\left({\xi }_3\right)}$>

进行积分；

步骤九：在几分过程中，如果不满足式

max(|h_lmn(t)|)≥|h₁₀₀(t)|/θ₁    l≠1  or

max(|v_lmn(t)|)≥|v₀₁₀(t)|/θ₁    m≠1  or，则继续积分；如果满足，则转步骤五；

max(|γ_lmn(t)|)≥|γ₀₀₁(t)|/θ₁  n≠1

以此类推，直至结束条件满足，即到300s结束。

步骤十：根据数学期望和方差的定义，求解状态均值和方差。

至此，本实例完毕。