CN109583007B - 一种火星进入飞行状态不确定性量化方法 - Google Patents

一种火星进入飞行状态不确定性量化方法 Download PDF

Info

Publication number
CN109583007B
CN109583007B CN201811212413.0A CN201811212413A CN109583007B CN 109583007 B CN109583007 B CN 109583007B CN 201811212413 A CN201811212413 A CN 201811212413A CN 109583007 B CN109583007 B CN 109583007B
Authority
CN
China
Prior art keywords
mars
uncertainty
state
random
polynomial
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811212413.0A
Other languages
English (en)
Other versions
CN109583007A (zh
Inventor
李爽
江秀强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN201811212413.0A priority Critical patent/CN109583007B/zh
Publication of CN109583007A publication Critical patent/CN109583007A/zh
Application granted granted Critical
Publication of CN109583007B publication Critical patent/CN109583007B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Operations Research (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种火星进入飞行状态不确定性量化方法，通过多项式混沌理论将随机火星进入动力学转化为等效的确定性动力学微分方程组。利用谱分解理论实现在火星进入飞行状态不确定性的统计特征剧烈变化时对多项式混沌中的正交基的更新，从而保证在多项式阶数不增加的情况下对不确定性量化的准确度。采用随机空间分解理论抑制在初始状态不确定性情况下的方程维数的剧烈增加，从而提高计算效率。本发明适用于对火星进入飞行的瞬时状态不确定性量化和状态演化不确定性量化。

Description

一种火星进入飞行状态不确定性量化方法
技术领域
本发明属于飞行力学技术领域，具体涉及一种火星进入飞行状态不确定性量化方法。
背景技术
火星进入段起始于飞行器到达火星大气层边界，终止于降落伞顺利展开，是整个火星探测任务中最关键的阶段之一。火星进入飞行技术是火星着陆探测任务的关键技术之一。然而，存在于火星大气密度、飞行器气动参数和初始状态中的不确定性通常会导致较大的飞行轨迹偏差和着陆偏差。因此，对火星进入状态轨迹的不确定性进行量化评估，以便于理清火星进入不确定性演化规律，是开展火星大气进入飞行系统设计和研究的重要前提。例如，鲁棒优化设计、风险与可靠性评估、以及任务前的高保真仿真与系统验证等都需要提前知道火星进入不确定性及其演化规律，避免单纯依赖经验进行相关设计。但现有的不确定性量化方法大都集中于流体力学、结构力学、多体动力学等领域，尚缺乏适用于火星进入飞行状态不确定性量化的方法。并且由于火星进入动力学具有高维和强非线性特征，难以直接应用传统的不确定性量化方法实现对火星进入不确定性的快速准确地量化评估。
不确定性量化是定量描述和评估系统中不确定性及其演化的一种手段。蒙特卡洛仿真是一种可用于分析不确定性分布的常用技术，它利用大量的随机样本作为初始输入条件，按照给定的步长进行演化计算，得出动力学系统响应，进而分析出状态轨迹的统计特性。但由于火星进入动力学是一个高维非线性系统，且不确定性参数较多，蒙特卡洛仿真的效率低下，为了得到较准确的分析结果需要大量的仿真计算次数。线性协方差分析法也常被用于线性或弱非线性系统的不确定性演化分析，但由于火星进入动力学系统具有全局的高维强非线性，线性协方差分析所基于的局部线性化难以准确近似原系统，这将导致不确定性量化分析难以达到较好的精度。多项式混沌是近年来研究不确定性量化的一类热点方法，其基本思想是利用正交多项式处理带有随机变量的随机过程问题，将系统的解在随机参数空间进行有限阶的多项式展开，通过伽辽金投影得出关于展开系数的方程组，进而求解出包含统计特性的精确解，从而以显式表达不确定性随状态轨迹的演化。但由于传统的多项式混沌需要对原系统模型进行改写，这对于具有高维状态和高维不确定性参数的火星进入不确定性量化问题而言，联立的方程组数量将远大于原问题的维数，使得求解效率大大降低。
发明内容
针对于上述现有技术的不足，本发明的目的在于提供一种火星进入飞行状态不确定性量化方法，以解决现有技术中火星进入动力学具有高维和强非线性特征，难以直接应用传统的不确定性量化方法实现对火星进入飞行状态不确定性的快速准确地量化评估的问题。
本发明方法的原理是利用随机变量表征火星进入动力学的不确定性，结合多项式混沌计算框架，通过考虑不确定性的统计特征随动力学非线性的变化来提高传统的不确定性量化方法的精度，通过对随机空间的自适应分解来抑制等效方程组维数的增加，从而得到一种适用于高维非线性火星进入动力学系统的不确定性量化方法，以实现高效准确的火星进入飞行状态不确定性量化。
为达到上述目的，本发明采用的技术方案如下：
本发明的一种火星进入飞行状态不确定性量化方法，包括步骤如下：
1)建立火星进入动力学模型；
2)建立含有不确定性的火星大气进入随机动力学模型；
3)获得等效的确定性微分方程组；
4)对等效的确定性微分方程组执行谱分解；
5)对系统状态执行随机空间分解；
6)获得火星进入飞行状态的均值和方差。
优选地，所述步骤1)具体包括：在火星质心固联坐标系下，考虑重力和气动力，则火星进入段飞行器的三自由度非线性动力学模型为：
Figure BDA0001832719660000026
Figure BDA0001832719660000021
Figure BDA0001832719660000022
其中，r为从火星质心到飞行器质心的距离，v为飞行器的速度，γ为飞行路径角，
Figure BDA0001832719660000023
为倾侧角，gm为火星重力加速度，L为飞行器受到的气动升力加速度，D为飞行器受到的气动阻力加速度；
Figure BDA0001832719660000024
式中，μ为火星引力常数，ρ为当地火星大气密度，S为飞行器的参考面积，m为飞行器质量，CL为飞行器的升力系数，CD为飞行器的阻力系数；弹道系数定义为：
Bc＝CDS/m (5)
升阻比定义为：
k＝CL/CD (6)。
优选地，所述步骤2)具体包括：引入随机变量表征火星进入动力学的不确定性，即：
a.初始状态不确定性模型为：
Figure BDA0001832719660000025
Figure BDA0001832719660000031
Figure BDA0001832719660000032
其中，随机变量ζ∈[-1,1]，ζ服从正态分布或均匀分布(即ζ的概率密度函数fζ(ζ)由用户根据ζ的分布类型来设定)；
Figure BDA0001832719660000033
为初始高度的标称值，/>
Figure BDA0001832719660000034
为初始速度的标称值，/>
Figure BDA0001832719660000035
为初始飞行路径角的标称值；/>
Figure BDA00018327196600000316
为初始高度的标准差，/>
Figure BDA00018327196600000317
为初始速度的标准差，/>
Figure BDA00018327196600000318
为初始飞行路径角的标准差；t0表示初始时刻；
b.动力学参数不确定性模型为：
Figure BDA0001832719660000036
Figure BDA0001832719660000037
Figure BDA0001832719660000038
其中，
Figure BDA0001832719660000039
为弹道系数的标称值，/>
Figure BDA00018327196600000310
为升阻比的标称值，/>
Figure BDA00018327196600000311
为火星大气密度比的标称值；/>
Figure BDA00018327196600000312
为弹道系数的标准差，δk为升阻比的标准差，δε为火星大气密度比的标准差；火星大气密度比为含有不确定偏差的火星大气密度与标称火星大气密度之比；
c.考虑动力学参数不确定性，则上述步骤1)中的火星进入动力学模型被改造成火星进入随机非线性动力学，即：
Figure BDA00018327196600000313
Figure BDA00018327196600000314
Figure BDA00018327196600000315
其中，h为飞行高度，即h＝r-r0，r0为当地的火星半径，ε为火星大气密度比。
优选地，所述步骤3)具体包括：
d.引入正交多项式Hm(ζ)作为基函数，并把步骤2)中的初始状态不确定性模型推广到任意时刻的状态形式：
Figure BDA0001832719660000041
Figure BDA0001832719660000042
Figure BDA0001832719660000043
其中，当ζ服从正态分布时，Hm(ζ)为Heimite多项式；当ζ服从均匀分布时，Hm(ζ)为Legendre多项式；hm(t),vm(t),γm(t)分别为h(t,ζ),v(t,ζ),γ(t,ζ)的多项式系数；t为时间变量；下标m为在闭区间[1,P]上的整数，P为用户自定义的多项式的最高阶数；
e.引入正交多项式H0(ζ)和H1(ζ)作为基函数，则步骤2)中的动力学参数不确定性模型被改写成：
Bc(ζ)＝Bc,0H0(ζ)+Bc,1H1(ζ) (19)
k(ζ)＝k0H0(ζ)+k1H1(ζ) (20)
ε(ζ)＝ε0H0(ζ)+ε1H1(ζ) (21)
其中，当ζ服从正态分布时，H0(ζ)和H1(ζ)为Heimite多项式；当ζ服从均匀分布时，H0(ζ)和H1(ζ)为Legendre多项式；Bc,0和Bc,1为Bc(ζ)分别在H0(ζ)和H1(ζ)基上的多项式系数；k0和k1为k(ζ)分别在H0(ζ)和H1(ζ)基上的多项式系数；ε0和ε1为ε(ζ)分别在H0(ζ)和H1(ζ)基上的多项式系数；
f.对步骤2)中的火星进入随机非线性动力学模型取在基函数Hm(ζ)上的伽辽金映射(Galerkin Projection)，得到在基函数Hm(ζ)上的等效确定性微分方程组：
Figure BDA0001832719660000044
Figure BDA0001832719660000045
Figure BDA0001832719660000051
其中，下标i,j,k都是在闭区间[0,P]上的自然数，且m,i,j,k的取值相互独立，<·>表示括号内表达式的数学期望；
Figure BDA0001832719660000052
则等效的确定性微分方程组(即公式(22)(23)(24))的解用来表征带有不确定性的火星进入飞行状态(即h(t,ζ)、v(t,ζ)、γ(t,ζ))；通过获得带有不确定性的火星进入飞行状态的统计特征参数，即可量化表示火星进入飞行状态不确定性。
优选地，所述步骤4)具体包括：
引入一组自适应判据，适时对等效确定性微分方程组执行谱分解；
g.自适应判据定义为：
Figure BDA0001832719660000053
Figure BDA0001832719660000054
Figure BDA0001832719660000055
其中，κ1,κ2,κ3为阈值实数(由用户根据状态误差发散程度来设定)，若自适应判据(即公式(26)(27)(28))中的任意一个成立，则执行等效确定性微分方程组的谱分解如下步骤h，若该自适应判据均不成立，则跳转至步骤5)；
h.引入三个新的随机变量ξ1,ξ2,ξ3作为新的正交基变量，即：
Figure BDA0001832719660000056
Figure BDA0001832719660000057
/>
Figure BDA0001832719660000058
其中，Z1,Z2,Z3分别表示ζ到ξ1,ξ2,ξ3的映射，P为用户自定义的多项式的最高阶数；ξj的概率密度函数
Figure BDA0001832719660000061
的计算公式为：
Figure BDA0001832719660000062
其中，下标j＝1,2,3；ζn为方程ξj-Zj(ζ)＝0的第n个实根，n为正整数；把ζn代入步骤2)中所设定的ζ的概率密度函数fζ(ζ)中，即得ζn的概率密度函数fζ(ζn)；
i.以ξ1,ξ2,ξ3为基变量，以概率密度函数
Figure BDA0001832719660000063
为权重系数，构建新的正交多项式/>
Figure BDA0001832719660000064
根据施密特正交化原则，验证/>
Figure BDA0001832719660000065
具有以下特性：
Figure BDA0001832719660000066
Figure BDA0001832719660000067
其中，Ci为任意实常数，δik为Kronecker-delta函数，l,m,n为在闭区间[1,P]上的整数，P为用户自定义的多项式的最高阶数；根据公式(33)和(34)求得新的正交多项式
Figure BDA0001832719660000068
j.步骤3)中公式(16)(17)(18)被更新为：
Figure BDA0001832719660000069
Figure BDA00018327196600000610
Figure BDA00018327196600000611
其中，
Figure BDA00018327196600000612
Figure BDA0001832719660000071
Figure BDA0001832719660000072
/>
于是步骤3)中的等效确定性微分方程组被更新为：
Figure BDA0001832719660000073
Figure BDA0001832719660000074
Figure BDA0001832719660000075
优选地，所述步骤5)具体包括：引入一组自适应判据用以决定是否执行随机空间分解；
k.记系统状态向量为x＝[h,v,γ]T，为便于以通式形式表述，用上标(ι)表示第ι个状态变量，ι＝1,2,3；则自适应判据被定义为：
Figure BDA0001832719660000076
Figure BDA0001832719660000077
其中，
Figure BDA0001832719660000078
α、θ1、θ2为自定义阈值参数，取值范围分别为α∈(0,1)，θ1∈(0,0.1)，θ2∈(0,0.1)；si表示第i个随机维度的灵敏度，sj表示第j个随机维度的灵敏度；下标i,p表示由p阶多项式表示的第i个随机维度，且p＝2,3,…,10；当i＝1,2,3时，φi分别代表步骤4)中的/>
Figure BDA0001832719660000079
Figure BDA00018327196600000710
表示进行随机空间分解之前的状态变量；Bk＝[(ai,bi)]3×1,i＝1,2,3，ai,bi分别表示第i个初始状态变量的上界和下界；若该自适应判据不成立，则跳转至步骤6)；若该自适应判据成立，则执行如下计算：
l.随机空间分解；若自适应判据被满足，随机域上的状态变量由新的基表达为：
Figure BDA0001832719660000081
更新的随机域写成：
Figure BDA0001832719660000082
其中，随机变量
Figure BDA0001832719660000083
为了确定系数/>
Figure BDA0001832719660000084
在[-1,1]3上选择(Np+1)个网格点/>
Figure BDA0001832719660000085
Np为[1,10]上的一个整数，i为整数且i∈[0,Np]；通过求解状态空间方程即下述公式(48)得系数/>
Figure BDA0001832719660000086
/>
Figure BDA0001832719660000087
其中
Figure BDA0001832719660000088
j为整数且j∈[0,Np]。
优选地，所述步骤6)具体包括：根据步骤5)的计算结果，火星进入飞行第ι个状态均值的计算公式为：
Figure BDA0001832719660000089
其中，ι＝1,2,3；
火星进入飞行第ι个状态的方差的计算公式为：
Figure BDA00018327196600000810
其中，
Figure BDA00018327196600000811
本发明的有益效果：
(1)本发明专门针对火星进入飞行状态不确定性量化问题，同时考虑了火星进入的初始状态不确定性和火星动力学参数不确定性，以及火星进入动力学的高维非线性特征，使得所得结果更符合实际。其中，火星进入动力学参数不确定性包含了火星大气密度不确定性、飞行器弹道系数不确定性、飞行器升阻比不确定性。
(2)本发明充分发挥多项式混沌、谱分解和随机空间分解方法的综合性能优势，在保证对火星进入飞行状态不确定性量化的较高精度的同时也提高了计算效率。
(3)本发明同时适用于对火星进入飞行的瞬时状态不确定性量化和状态演化不确定性量化。
附图说明
图1为本发明方法的原理流程图。
图2为火星进入的高度与速度不确定性量化结果示意图。
图3为火星进入的高度与飞行路径角不确定性量化结果示意图。
图4为火星进入的速度与飞行路径角不确定性量化结果示意图。
图5为火星进入的高度方差对比示意图。
图6为火星进入的速度方差对比示意图。
图7为火星进入的飞行路径角方差对比示意图。
图8为火星进入的高度均值误差对比示意图。
图9为火星进入的速度均值误差对比示意图。
图10为火星进入的飞行路径角均值误差对比示意图。
具体实施方式
为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。
本发明首先建立火星初始状态不确定性和动力学参数不确定性的随机数学模型，设置随机数学模型中各项的概率分布，使得传统的确定性的火星进入动力学被改造成随机的火星进入动力学。然后，利用多项式混沌展开式将随机的火星进入动力学转化为以正交多项式为基的等效的确定性微分方程组。接着，根据火星进入状态不确定性受火星进入动力学非线性的影响而变化的剧烈程度，利用谱分解方法对正交多项式基进行更新，以保证由多项式混沌转化的等效确定性微分方程组的解具有较好的精度。同时，对输入的初始状态不确定性进行随机空间分解，以抑制等效确定性微分方程组的维数的增加，提高计算效率。等效确定性微分方程组的解就是含有不确定性的火星进入飞行状态。最后，根据含有不确定性的火星进入飞行状态可以得到火星进入飞行状态的均值和方差，进而最终实现对火星进入飞行状态的不确定性量化。具体步骤如下：
参照图1所示，本发明的一种火星进入飞行状态不确定性量化方法，包括步骤如下：
1)建立火星进入动力学模型；
在火星质心固联坐标系下，考虑重力和气动力，则火星进入段飞行器的三自由度非线性动力学模型为：
Figure BDA0001832719660000101
Figure BDA0001832719660000102
Figure BDA0001832719660000103
其中，r为从火星质心到飞行器质心的距离，v为飞行器的速度，γ为飞行路径角，
Figure BDA0001832719660000104
为倾侧角，gm为火星重力加速度，L为飞行器受到的气动升力加速度，D为飞行器受到的气动阻力加速度；
Figure BDA0001832719660000105
式中，μ为火星引力常数，ρ为当地火星大气密度，S为飞行器的参考面积，m为飞行器质量，CL为飞行器的升力系数，CD为飞行器的阻力系数；弹道系数定义为：
Bc＝CDS/m (5)
升阻比定义为：
k＝CL/CD (6)。
2)建立含有不确定性的火星大气进入随机动力学模型；
引入随机变量表征火星进入动力学的不确定性，即：
a.初始状态不确定性模型为：
Figure BDA0001832719660000106
Figure BDA0001832719660000107
Figure BDA0001832719660000108
其中，随机变量ζ∈[-1,1]，ζ服从正态分布或均匀分布(即ζ的概率密度函数fζ(ζ)由用户根据ζ的分布类型来设定)；
Figure BDA0001832719660000109
为初始高度的标称值，/>
Figure BDA00018327196600001010
为初始速度的标称值，/>
Figure BDA00018327196600001011
为初始飞行路径角的标称值；/>
Figure BDA00018327196600001012
为初始高度的标准差，/>
Figure BDA00018327196600001013
为初始速度的标准差，/>
Figure BDA00018327196600001014
为初始飞行路径角的标准差；t0表示初始时刻；
b.动力学参数不确定性模型为：
Figure BDA0001832719660000111
Figure BDA0001832719660000112
Figure BDA0001832719660000113
其中，
Figure BDA0001832719660000114
为弹道系数的标称值，/>
Figure BDA0001832719660000115
为升阻比的标称值，/>
Figure BDA0001832719660000116
为火星大气密度比的标称值；/>
Figure BDA0001832719660000117
为弹道系数的标准差，δk为升阻比的标准差，δε为火星大气密度比的标准差；火星大气密度比为含有不确定偏差的火星大气密度与标称火星大气密度之比；
c.考虑动力学参数不确定性，则上述步骤1)中的火星进入动力学模型被改造成火星进入随机非线性动力学，即：
Figure BDA0001832719660000118
Figure BDA0001832719660000119
Figure BDA00018327196600001110
其中，h为飞行高度，即h＝r-r0，r0为当地的火星半径，ε为火星大气密度比。
3)获得等效的确定性微分方程组；
d.引入正交多项式Hm(ζ)作为基函数，并把步骤2)中的初始状态不确定性模型推广到任意时刻的状态形式：
Figure BDA00018327196600001111
Figure BDA00018327196600001112
Figure BDA00018327196600001113
其中，当ζ服从正态分布时，Hm(ζ)为Heimite多项式；当ζ服从均匀分布时，Hm(ζ)为Legendre多项式；hm(t),vm(t),γm(t)分别为h(t,ζ),v(t,ζ),γ(t,ζ)的多项式系数；t为时间变量；下标m为在闭区间[1,P]上的整数，P为用户自定义的多项式的最高阶数；
e.引入正交多项式H0(ζ)和H1(ζ)作为基函数，则步骤2)中的动力学参数不确定性模型被改写成：
Bc(ζ)＝Bc,0H0(ζ)+Bc,1H1(ζ) (19)
k(ζ)＝k0H0(ζ)+k1H1(ζ) (20)
ε(ζ)＝ε0H0(ζ)+ε1H1(ζ) (21)
其中，当ζ服从正态分布时，H0(ζ)和H1(ζ)为Heimite多项式；当ζ服从均匀分布时，H0(ζ)和H1(ζ)为Legendre多项式；Bc,0和Bc,1为Bc(ζ)分别在H0(ζ)和H1(ζ)基上的多项式系数；k0和k1为k(ζ)分别在H0(ζ)和H1(ζ)基上的多项式系数；ε0和ε1为ε(ζ)分别在H0(ζ)和H1(ζ)基上的多项式系数；
f.对步骤2)中的火星进入随机非线性动力学模型取在基函数Hm(ζ)上的伽辽金映射(Galerkin Projection)，得到在基函数Hm(ζ)上的等效确定性微分方程组：
Figure BDA0001832719660000121
Figure BDA0001832719660000122
Figure BDA0001832719660000123
其中，下标i,j,k都是在闭区间[0,P]上的自然数，且m,i,j,k的取值相互独立，<·>表示括号内表达式的数学期望；
Figure BDA0001832719660000131
则等效的确定性微分方程组(即公式(22)(23)(24))的解就可以用来表征带有不确定性的火星进入飞行状态(即h(t,ζ)、v(t,ζ)、γ(t,ζ))。通过获得带有不确定性的火星进入飞行状态的统计特征参数(即均值和方差)，即可量化表示火星进入飞行状态不确定性。
4)对等效确定性微分方程组执行谱分解；
为了克服由飞行状态不确定性的统计特征的急剧变化引起的发散问题，引入一组自适应判据，适时对等效确定性微分方程组执行谱分解；
g.自适应判据定义为：
Figure BDA0001832719660000132
Figure BDA0001832719660000133
/>
Figure BDA0001832719660000134
其中，κ1,κ2,κ3为阈值实数(由用户根据状态误差发散程度来设定)。若自适应判据(即公式(26)(27)(28))中的任意一个成立，则执行等效确定性微分方程组的谱分解如下步骤h，若该自适应判据均不成立，则跳转至步骤5)；
h.引入三个新的随机变量ξ1,ξ2,ξ3作为新的正交基变量，即：
Figure BDA0001832719660000135
Figure BDA0001832719660000136
Figure BDA0001832719660000137
其中，Z1,Z2,Z3分别表示ζ到ξ1,ξ2,ξ3的映射，P为用户自定义的多项式的最高阶数；ξj的概率密度函数
Figure BDA0001832719660000138
的计算公式为：
Figure BDA0001832719660000139
其中，下标j＝1,2,3；ζn为方程ξj-Zj(ζ)＝0的第n个实根，n为正整数；把ζn代入步骤2)中所设定的ζ的概率密度函数fζ(ζ)中，即得ζn的概率密度函数fζ(ζn)；
i.以ξ1,ξ2,ξ3为基变量，以概率密度函数
Figure BDA0001832719660000141
为权重系数，构建新的正交多项式/>
Figure BDA0001832719660000142
根据施密特正交化原则，验证/>
Figure BDA0001832719660000143
具有以下特性：
Figure BDA0001832719660000144
Figure BDA0001832719660000145
其中，Ci为任意实常数，δik为Kronecker-delta函数，l,m,n为在闭区间[1,P]上的整数，P为用户自定义的多项式的最高阶数；根据公式(33)和(34)求得新的正交多项式
Figure BDA0001832719660000146
j.步骤3)中公式(16)(17)(18)被更新为：
Figure BDA0001832719660000147
Figure BDA0001832719660000148
Figure BDA0001832719660000149
其中，
Figure BDA00018327196600001410
Figure BDA00018327196600001411
Figure BDA00018327196600001412
于是步骤3)中的等效确定性微分方程组被更新为：
Figure BDA0001832719660000151
Figure BDA0001832719660000152
Figure BDA0001832719660000153
5)对系统状态执行随机空间分解；
为了抑制初始状态不确定性量化精度受动力学参数不确定性和系统非线性的影响而降低，同时抑制计算量的剧烈增涨，适时对初始状态所在的随机空间进行分解，并引入一组自适应判据用以决定是否执行随机空间分解；
k.记系统状态向量为x＝[h,v,γ]T，为便于以通式形式表述，用上标(ι)表示第ι个状态变量，ι＝1,2,3；则自适应判据被定义为：
Figure BDA0001832719660000154
Figure BDA0001832719660000155
其中，
Figure BDA0001832719660000156
α、θ1、θ2为自定义阈值参数，取值范围分别为α∈(0,1)，θ1∈(0,0.1)，θ2∈(0,0.1)；si表示第i个随机维度的灵敏度，sj表示第j个随机维度的灵敏度；下标i,p表示由p阶多项式表示的第i个随机维度，且p＝2,3,…,10；当i＝1,2,3时，φi分别代表步骤4)中的/>
Figure BDA0001832719660000157
Figure BDA0001832719660000158
表示进行随机空间分解之前的状态变量；Bk＝[(ai,bi)]3×1,i＝1,2,3，ai,bi分别表示第i个初始状态变量的上界和下界；若该自适应判据不成立，则跳转至步骤6)；若该自适应判据成立，则执行如下计算：
l.随机空间分解；若自适应判据被满足，随机域上的状态变量由新的基表达为：
Figure BDA0001832719660000159
更新的随机域写成：
Figure BDA0001832719660000161
其中，随机变量
Figure BDA0001832719660000162
为了确定系数/>
Figure BDA0001832719660000163
在[-1,1]3上选择(Np+1)个网格点/>
Figure BDA0001832719660000164
Np为[1,10]上的一个整数，i为整数且i∈[0,Np]；通过求解状态空间方程即下述公式(48)得系数/>
Figure BDA0001832719660000165
Figure BDA0001832719660000166
其中
Figure BDA0001832719660000167
j为整数且j∈[0,Np]。
6)获得火星进入飞行状态的均值和方差；
根据步骤5)的计算结果，火星进入飞行第ι个状态均值的计算公式为：
Figure BDA0001832719660000168
其中，ι＝1,2,3；
火星进入飞行第ι个状态的方差的计算公式为：
Figure BDA0001832719660000169
其中，
Figure BDA00018327196600001610
本发明方法的实例：结合图2和图3说明本发明的实例验证，设定如下计算条件和技术参数。即，根据给定的初始状态和动力学参数的概率分布及其参数，在速度v＞450m/s且高度h＞6.5km的飞行阶段，对火星进入飞行状态轨迹进行不确定性量化。
公式(7)(8)(9)(10)(11)(12)所表达的初始状态和动力学参数不确定性的概率分布分别设置为如表1所示：
表1
参数 均值 高斯分布3σ方差 单位
初始高度,h0 125 12.5 km
初始速度,v0 5500 550 m/s
初始飞行路径角,γ0 -15.2 1.52 deg
弹道系数,Bc 121.6 12.16 kg/m2
升阻比,k 0.24 0.024 -
大气密度不确定性,ε 1 0.1 -
公式(26)(27)(28)(44)(45)所表达的自适应判据中的参数分布设置为：κ1＝κ2＝κ3＝2/3，θ1＝10-2，θ2＝10-1。步骤3和步骤4中的多项式最高阶数P取4。
实例所运行的软件环境为
Figure BDA0001832719660000171
(2016b,64位版)，硬件条件为Intel Core(TM)i7-7820HQ CPU@2.9GHz和32GB内存的/>
Figure BDA0001832719660000172
移动工作站。
图2为火星进入的高度与速度不确定性量化结果。图3为火星进入的高度与飞行路径角不确定性量化结果。图4为火星进入的速度与飞行路径角不确定性量化结果。为了更好的说明本发明方法的优势，采用传统的蒙特卡洛(Monte-Carlo)方法和一般的多项式混沌方法的运行结果与本发明方法的运行结果进行对比，对比结果如图5、图6、图7、图8、图9、图10所示。图5为火星进入的高度方差对比。图6为火星进入的速度方差对比。图7为火星进入的飞行路径角方差对比。图8为火星进入的高度均值误差对比。图9为火星进入的速度均值误差对比。图10为火星进入的飞行路径角均值误差对比。从图2、图3、图4可以看出，本发明方法能够实现对火星进入飞行状态的不确定性量化。从图5、图6、图7、图8、图9、图10和表2所示的对比结果可以看出：(1)为获得同等精度量级的不确定性量化结果，本发明方法的运行时间不到蒙特卡洛方法的运行时间的5％；(2)一般的多项式混沌方法比本发明多花30多秒的时间，但本发明方法的结果误差比一般的多项式混沌方法的结果误差小1～2个数量级。结果综合表明，本发明方法能够高效进行火星进入飞行状态的不确定性量化。表2如下：
表2
蒙特卡洛方法 一般的多项式混沌方法 本发明方法 单位
3595.849 203.695 171.925 秒
本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。

Claims (5)

1.一种火星进入飞行状态不确定性量化方法，其特征在于，包括步骤如下：
1)建立火星进入动力学模型；
2)建立含有不确定性的火星大气进入随机动力学模型；
3)获得等效的确定性微分方程组；
4)对等效的确定性微分方程组执行谱分解；
5)对系统状态执行随机空间分解；
6)获得火星进入飞行状态的均值和方差；
所述步骤1)具体包括：在火星质心固联坐标系下，考虑重力和气动力，则火星进入段飞行器的三自由度非线性动力学模型为：
Figure FDA0003917590330000011
Figure FDA0003917590330000012
Figure FDA0003917590330000013
其中，r为从火星质心到飞行器质心的距离，v为飞行器的速度，γ为飞行路径角，
Figure FDA0003917590330000014
为倾侧角，gm为火星重力加速度，L为飞行器受到的气动升力加速度，D为飞行器受到的气动阻力加速度；
gm＝μ/r2，
Figure FDA0003917590330000015
式中，μ为火星引力常数，ρ为当地火星大气密度，S为飞行器的参考面积，m为飞行器质量，CL为飞行器的升力系数，CD为飞行器的阻力系数；弹道系数定义为：
Bc＝CDS/m (5)
升阻比定义为：
k＝CL/CD (6)；
所述步骤2)具体包括：引入随机变量表征火星进入动力学的不确定性，即：
a.初始状态不确定性模型为：
Figure FDA0003917590330000016
Figure FDA0003917590330000017
Figure FDA0003917590330000018
其中，随机变量ζ∈[-1,1]，ζ服从正态分布或均匀分布；
Figure FDA0003917590330000019
为初始高度的标称值，/>
Figure FDA00039175903300000110
为初始速度的标称值，/>
Figure FDA0003917590330000021
为初始飞行路径角的标称值；/>
Figure FDA0003917590330000022
为初始高度的标准差，/>
Figure FDA0003917590330000023
为初始速度的标准差，/>
Figure FDA0003917590330000024
为初始飞行路径角的标准差；t0表示初始时刻；
b.动力学参数不确定性模型为：
Figure FDA0003917590330000025
Figure FDA0003917590330000026
Figure FDA0003917590330000027
其中，
Figure FDA0003917590330000028
为弹道系数的标称值，/>
Figure FDA0003917590330000029
为升阻比的标称值，/>
Figure FDA00039175903300000210
为火星大气密度比的标称值；/>
Figure FDA00039175903300000211
为弹道系数的标准差，δk为升阻比的标准差，δε为火星大气密度比的标准差；火星大气密度比为含有不确定偏差的火星大气密度与标称火星大气密度之比；
c.考虑动力学参数不确定性，则上述步骤1)中的火星进入动力学模型被改造成火星进入随机非线性动力学，即：
Figure FDA00039175903300000212
Figure FDA00039175903300000213
Figure FDA00039175903300000214
其中，h为飞行高度，即h＝r-r0，r0为当地的火星半径，ε为火星大气密度比。
2.根据权利要求1所述的火星进入飞行状态不确定性量化方法，其特征在于，所述步骤3)具体包括：
d.引入正交多项式Hm(ζ)作为基函数，并把步骤2)中的初始状态不确定性模型推广到任意时刻的状态形式：
Figure FDA00039175903300000215
Figure FDA00039175903300000216
Figure FDA0003917590330000031
其中，当ζ服从正态分布时，Hm(ζ)为Heimite多项式；当ζ服从均匀分布时，Hm(ζ)为Legendre多项式；hm(t),vm(t),γm(t)分别为h(t,ζ),v(t,ζ),γ(t,ζ)的多项式系数；t为时间变量；下标m为在闭区间[1,P]上的整数，P为用户自定义的多项式的最高阶数；
e.引入正交多项式H0(ζ)和H1(ζ)作为基函数，则步骤2)中的动力学参数不确定性模型被改写成：
Bc(ζ)＝Bc,0H0(ζ)+Bc,1H1(ζ) (19)
k(ζ)＝k0H0(ζ)+k1H1(ζ) (20)
ε(ζ)＝ε0H0(ζ)+ε1H1(ζ) (21)
其中，当ζ服从正态分布时，H0(ζ)和H1(ζ)为Heimite多项式；当ζ服从均匀分布时，H0(ζ)和H1(ζ)为Legendre多项式；Bc,0和Bc,1为Bc(ζ)分别在H0(ζ)和H1(ζ)基上的多项式系数；k0和k1为k(ζ)分别在H0(ζ)和H1(ζ)基上的多项式系数；ε0和ε1为ε(ζ)分别在H0(ζ)和H1(ζ)基上的多项式系数；
f.对步骤2)中的火星进入随机非线性动力学模型取在基函数Hm(ζ)上的伽辽金映射，得到在基函数Hm(ζ)上的等效确定性微分方程组：
Figure FDA0003917590330000032
/>
Figure FDA0003917590330000033
Figure FDA0003917590330000034
其中，下标i,j,k都是在闭区间[0,P]上的自然数，且m,i,j,k的取值相互独立，<·>表示括号内表达式的数学期望；
Figure FDA0003917590330000041
则等效的确定性微分方程组的解用来表征带有不确定性的火星进入飞行状态；通过获得带有不确定性的火星进入飞行状态的统计特征参数，即可量化表示火星进入飞行状态不确定性。
3.根据权利要求2所述的火星进入飞行状态不确定性量化方法，其特征在于，所述步骤4)具体包括：
引入一组自适应判据，适时对等效确定性微分方程组执行谱分解；
g.自适应判据定义为：
Figure FDA0003917590330000042
Figure FDA0003917590330000043
Figure FDA0003917590330000044
其中，κ1,κ2,κ3为阈值实数，若自适应判据中的任意一个成立，则执行等效确定性微分方程组的谱分解如下步骤h，若该自适应判据均不成立，则跳转至步骤5)；
h.引入三个新的随机变量ξ1,ξ2,ξ3作为新的正交基变量，即：
Figure FDA0003917590330000045
Figure FDA0003917590330000046
Figure FDA0003917590330000047
其中，Z1,Z2,Z3分别表示ζ到ξ1,ξ2,ξ3的映射，P为用户自定义的多项式的最高阶数；ξj的概率密度函数
Figure FDA0003917590330000048
的计算公式为：/>
Figure FDA0003917590330000049
其中，下标j＝1,2,3；ζn为方程ξj-Zj(ζ)＝0的第n个实根，n为正整数；把ζn代入步骤2)中所设定的ζ的概率密度函数fζ(ζ)中，即得ζn的概率密度函数fζ(ζn)；
i.以ξ1,ξ2,ξ3为基变量，以概率密度函数
Figure FDA0003917590330000051
为权重系数，构建新的正交多项式/>
Figure FDA0003917590330000052
根据施密特正交化原则，验证/>
Figure FDA0003917590330000053
具有以下特性：
Figure FDA0003917590330000054
Figure FDA0003917590330000055
其中，Ci为任意实常数，δik为Kronecker-delta函数，l,m,n为在闭区间[1,P]上的整数，P为用户自定义的多项式的最高阶数；根据公式(33)和(34)求得新的正交多项式
Figure FDA0003917590330000056
j.步骤3)中公式(16)(17)(18)被更新为：
Figure FDA0003917590330000057
Figure FDA0003917590330000058
Figure FDA0003917590330000059
其中，
Figure FDA00039175903300000510
Figure FDA00039175903300000511
Figure FDA00039175903300000512
于是步骤3)中的等效确定性微分方程组被更新为：
Figure FDA00039175903300000513
Figure FDA0003917590330000061
/>
Figure FDA0003917590330000062
4.根据权利要求3所述的火星进入飞行状态不确定性量化方法，其特征在于，所述步骤5)具体包括：引入一组自适应判据用以决定是否执行随机空间分解；
k.记系统状态向量为x＝[h,v,γ]T，为便于以通式形式表述，用上标(ι)表示第ι个状态变量，ι＝1,2,3；则自适应判据被定义为：
Figure FDA0003917590330000063
Figure FDA0003917590330000064
其中，
Figure FDA0003917590330000065
α、θ1、θ2为自定义阈值参数，取值范围分别为α∈(0,1)，θ1∈(0,0.1)，θ2∈(0,0.1)；si表示第i个随机维度的灵敏度，sj表示第j个随机维度的灵敏度；下标i,p表示由p阶多项式表示的第i个随机维度，且p＝2,3,…,10；当i＝1,2,3时，φi分别代表步骤4)中的/>
Figure FDA0003917590330000066
Figure FDA00039175903300000612
表示进行随机空间分解之前的状态变量；Bk＝[(ai,bi)]3×1,i＝1,2,3，ai,bi分别表示第i个初始状态变量的上界和下界；若该自适应判据不成立，则跳转至步骤6)；若该自适应判据成立，则执行如下计算：
l.随机空间分解；若自适应判据被满足，随机域上的状态变量由新的基表达为：
Figure FDA0003917590330000067
更新的随机域写成：
Figure FDA0003917590330000068
其中，随机变量
Figure FDA0003917590330000069
为了确定系数/>
Figure FDA00039175903300000610
在[-1,1]3上选择(Np+1)个网格点/>
Figure FDA00039175903300000611
Np为[1,10]上的一个整数，i为整数且i∈[0,Np]；通过求解状态空间方程即下述公式(48)得系数/>
Figure FDA0003917590330000071
Figure FDA0003917590330000072
其中
Figure FDA0003917590330000073
j为整数且j∈[0,Np]。/>
5.根据权利要求4所述的火星进入飞行状态不确定性量化方法，其特征在于，所述步骤6)具体包括：根据步骤5)的计算结果，火星进入飞行第ι个状态均值的计算公式为：
Figure FDA0003917590330000074
其中，ι＝1,2,3；
火星进入飞行第ι个状态的方差的计算公式为：
Figure FDA0003917590330000075
其中，
Figure FDA0003917590330000076
/>
CN201811212413.0A 2018-10-12 2018-10-12 一种火星进入飞行状态不确定性量化方法 Active CN109583007B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811212413.0A CN109583007B (zh) 2018-10-12 2018-10-12 一种火星进入飞行状态不确定性量化方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811212413.0A CN109583007B (zh) 2018-10-12 2018-10-12 一种火星进入飞行状态不确定性量化方法

Publications (2)

Publication Number Publication Date
CN109583007A CN109583007A (zh) 2019-04-05
CN109583007B true CN109583007B (zh) 2023-05-26

Family

ID=65920508

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811212413.0A Active CN109583007B (zh) 2018-10-12 2018-10-12 一种火星进入飞行状态不确定性量化方法

Country Status (1)

Country Link
CN (1) CN109583007B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110442117B (zh) * 2019-08-26 2020-12-25 北京理工大学 一种火星探测器大底气动耦合分离过程安全性分析方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102930166A (zh) * 2012-11-05 2013-02-13 北京理工大学 基于混沌多项式的行星大气进入状态不确定度获取方法
CN106295000A (zh) * 2016-08-10 2017-01-04 北京理工大学 一种考虑不确定性影响的火星大气进入段轨迹优化方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102930166A (zh) * 2012-11-05 2013-02-13 北京理工大学 基于混沌多项式的行星大气进入状态不确定度获取方法
CN106295000A (zh) * 2016-08-10 2017-01-04 北京理工大学 一种考虑不确定性影响的火星大气进入段轨迹优化方法

Also Published As

Publication number Publication date
CN109583007A (zh) 2019-04-05

Similar Documents

Publication Publication Date Title
Yang et al. Some issues in uncertainty quantification and parameter tuning: A case study of convective parameterization scheme in the WRF regional climate model
Zheng et al. Inverse calculation approaches for source determination in hazardous chemical releases
da Silva et al. Ensemble-based state estimator for aerodynamic flows
CN108959182B (zh) 基于高斯过程回归的小天体引力场建模方法
CN112114521A (zh) 航天器智能预测控制进入制导方法
CN110059439B (zh) 一种基于数据驱动的航天器轨道确定方法
Metref et al. A non-Gaussian analysis scheme using rank histograms for ensemble data assimilation
CN114168796B (zh) 一种建立飞行器高空气动力数据库的方法
CN107402903A (zh) 基于微分代数与高斯和的非线性系统状态偏差演化方法
Abramov et al. A new algorithm for low-frequency climate response
CN111415010A (zh) 一种基于贝叶斯神经网络的风电机组参数辨识方法
CN104266650A (zh) 一种基于采样点继承策略的火星着陆器大气进入段导航方法
CN109583007B (zh) 一种火星进入飞行状态不确定性量化方法
CN112819303A (zh) 基于pce代理模型的飞行器追踪效能评估方法及系统
Chen et al. Proactive quality control: Observing system simulation experiments with the Lorenz’96 model
CN107766588B (zh) 逃逸飞行器遵循多种概率分布的多次碰撞情况仿真方法
CN113486591B (zh) 一种卷积神经网络结果的重力多参量数据密度加权反演方法
CN102930166B (zh) 基于混沌多项式的行星大气进入状态不确定度获取方法
CN105975677A (zh) 一种复杂外形低轨航天器气动特性的快速预测方法
Ray et al. Bayesian calibration of a RANS model with a complex response surface-a case study with jet-in-crossflow configuration
Cao et al. Flight Trajectory Simulation and Aerodynamic Parameter Identification of Large‐Scale Parachute
Jin et al. Monte Carlo simulation for aerodynamic coefficients of satellites in Low-Earth Orbit
Wang et al. Real-time data driven simulation of air contaminant dispersion using particle filter and UAV sensory system
CN116204985A (zh) 飞行器壁面应力张量抽样方法、装置、设备及存储介质
CN107609250B (zh) 地形重力波举力参数化方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant