CN107402903B - 基于微分代数与高斯和的非线性***状态偏差演化方法 - Google Patents

基于微分代数与高斯和的非线性***状态偏差演化方法 Download PDF

Info

Publication number
CN107402903B
CN107402903B CN201710551107.9A CN201710551107A CN107402903B CN 107402903 B CN107402903 B CN 107402903B CN 201710551107 A CN201710551107 A CN 201710551107A CN 107402903 B CN107402903 B CN 107402903B
Authority
CN
China
Prior art keywords
gaussian distribution
sub
deviation
matrix
gaussian
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
CN201710551107.9A
Other languages
English (en)
Other versions
CN107402903A (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.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
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 National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN201710551107.9A priority Critical patent/CN107402903B/zh
Publication of CN107402903A publication Critical patent/CN107402903A/zh
Application granted granted Critical
Publication of CN107402903B publication Critical patent/CN107402903B/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
    • 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
    • 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
    • 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/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Complex Calculations (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

本发明公开了一种基于微分代数与高斯和的非线性***状态偏差演化方法，步骤包括：采用微分代数方法预报非线性***的终端状态，并表示为关于初始状态偏差的高阶泰勒展开多项式；确定子高斯分布协方差矩阵，针对各个子高斯分布打靶拟合高斯和模型；计算子高斯分布高阶中心矩；确定各子高斯分布在终端时刻的均值和协方差矩阵，给出高斯和形式的终端状态偏差分布概率密度函数。本发明能够自动拓展至任意指定阶偏差演化精度，无需手动推导动力学方程的高阶偏导数，适用于非线性强和长时间预报的非线性***偏差演化分析问题，在精确描述偏差分布及其非高斯性的同时，仍能保持相对Monte Carlo仿真方法的显著效率优势，具有使用方便、计算精度高的优点。

Description

基于微分代数与高斯和的非线性***状态偏差演化方法
技术领域
本发明涉及航天器轨道动力学领域，具体涉及一种基于微分代数与高斯和的非线性***状态偏差演化方法，用于非线性***的状态及其偏差预报。
背景技术
航天器轨迹偏差演化技术是开展空间碰撞预警、航天器轨道确定、导航滤波以及航天器轨迹安全性分析与设计等航天任务的重要基础技术。轨迹偏差的预报精度直接影响了预警、定轨和导航结果的可信度，甚至决定着航天任务的成败。由于航天器轨道动力学模型的实质就是一个常微分方程组或状态转移方程形式的非线性***，因此航天器轨迹偏差演化问题实质就是非线性***状态偏差预报问题。
经典的非线性***偏差演化方法是线性化协方差分析方法和MonteCarlo仿真方法。线性化协方差分析法，先将非线性动力学***线性化，再采用线性***协方差分析方法进行偏差预报计算，该方法使用简单、计算量小，但是对强非线性***或长时间偏差演化预报问题计算误差较大。MonteCarlo仿真方法，通过统计多次打靶试验结果可以得到高精度的状态偏差演化结果，但是使用该方法需要进行大量抽样仿真，计算量大，很多情况下计算效率难以满足任务需求。
微分代数方法是一种在数值环境下计算非线性映射关于自变量的高阶偏导数的方法，利用这些高阶偏导数可以对非线性***进行任意阶次的泰勒展开，进而使用高阶多项式预报***状态变量。该方法最初由美国密歇根大学粒子物理学教授MartinBerz提出，基本原理可见参考文献：Berz M.Differential algebraic description of beamdynamics to very high orders[J].Part.Accel.,1988,24(SSC-152):109-124。目前国内尚未有应用该方法进行非线性动力学***偏差预报的研究或应用成果。
高斯和模型是以多个子高斯分布加权和的形式描述***状态偏差分布的方法，可用于准确描述非高斯分布的状态偏差概率密度函数。本发明面向一般的非线性动力学***，***动力学模型可以是状态转移方程形式，或常微分方程形式，利用微分代数方法求解终端状态关于初始状态偏差的高阶多项式，进而采用高斯和模型拟合初始状态偏差分布，并基于前述高阶多项式预报终端时刻的各子高斯分布，以此来准确描述终端状态的非高斯分布特性。
发明内容
本发明要解决的技术问题：针对现有技术的上述问题，提供一种能够自动拓展至任意指定阶偏差演化精度，无需手动推导动力学方程的高阶偏导数，适用于非线性强和长时间预报的非线性***偏差演化分析问题，在精确描述偏差分布及其非高斯性的同时，仍能保持相对MonteCarlo仿真方法的显著效率优势，使用方便、计算精度高的基于微分代数与高斯和的非线性***状态偏差演化方法。
为了解决上述技术问题，本发明采用的技术方案为：
一种基于微分代数与高斯和的非线性***状态偏差演化方法，用于航天器轨迹偏差演化，该方法的步骤包括：
1)输入非线性***的动力学方程、偏差演化阶数N、初始状态相对标称值
Figure GDA0002884026650000026
和初始偏差协方差矩阵P0；
2)采用微分代数方法，针对非线性***的动力学方程预报非线性***的终端状态，并表示为关于初始状态偏差δx0的高阶泰勒展开多项式；
3)针对初始偏差协方差矩阵P0确定子高斯分布协方差矩阵Pi，针对子高斯分布协方差矩阵Pi中的多个子高斯分布进行打靶拟合初始偏差分布生成高斯和模型；
4)针对生成高斯和模型计算给定阶数的子高斯分布高阶中心矩；
5)基于终端状态关于初始状态偏差的高阶多项式和子高斯分布的各阶中心矩，确定各子高斯分布在终端时刻的均值和协方差矩阵，进而给出高斯和形式的终端状态偏差分布概率密度函数。
优选地，步骤1)中非线性***的动力学方程为非线性映射形式xf＝f(x0)或常微分方程形式
Figure GDA0002884026650000021
且所述非线性映射形式的函数f(·)或者常微分方程形式的函数g(·)均为任意从n维实数域Rn到Rn的非线性映射。
优选地，步骤2)的详细步骤包括：
2.1)将初始状态x0初始化为式(1)所示包含初始状态相对标称值
Figure GDA0002884026650000022
的偏差形式；
Figure GDA0002884026650000023
式(1)中，x0为初始状态，
Figure GDA0002884026650000024
为初始状态相对标称值，δx0为初始状态偏差；
2.2)将初始状态x0代入非线性***的动力学方程，应用微分代数方法，求出关于初始状态偏差δx0的高阶泰勒展开多项式描述的终端状态如式(2)所示；
Figure GDA0002884026650000025
式(2)中，xf为终端状态，T为高阶泰勒展开多项式，T的上标N为偏差演化阶数、T的下标xf为高阶泰勒展开多项式描述的内容为终端状态，δx0为初始状态偏差。
优选地，步骤3)的详细步骤包括：
3.1)针对初始偏差协方差矩阵P0的逆进行Cholesky分解求平方根信息矩阵R0，且满足P0＝R0 -1R0 -T；
3.2)给定子高斯分布的平方根信息矩阵的下界Rmin，求子高斯分布协方差矩阵Pi；
3.3)将多个子高斯分布打靶生成高斯和模型来拟合初始高斯分布。
优选地，步骤3.2)中求子高斯分布协方差矩阵Pi的详细步骤包括：
3.2.1)采用式(3)所示函数表达式计算矩阵
Figure GDA0002884026650000031
的奇异值分解，其中R0为平方根信息矩阵，Rmin为子高斯分布的平方根信息矩阵的下界；
Figure GDA0002884026650000032
式(3)中，Ui和Vi为标准正交矩阵，对角阵Si＝diag(σi1,...,σin)，其中对角线元素σi1,...,σin是正奇异值，R0为平方根信息矩阵，Rmin为子高斯分布的平方根信息矩阵的下界；
3.2.2)对任意k＝1,...,n，判断对角阵Si中的第k个对角线元素σik大于或等于1是否成立，如果成立则所求子高斯分布的平方根信息矩阵Ri为平方根信息矩阵R0，跳转执行步骤3.2.5)；否则，跳转执行下一步；
3.2.3)建立式(4)所示n维对角矩阵，删除n维对角矩阵的全0行，得到奇异值变化矩阵δSi；
Figure GDA0002884026650000033
式(4)中，δSi_full为n维对角矩阵，σi1为对角阵Si中的第1个对角线元素，σin为对角阵Si中的第n个对角线元素；
3.2.4)根据式(5)求取信息变化矩阵δRi，并根据式(6)求解子高斯分布的平方根信息矩阵Ri；
δRi＝δSiVi TRmin (5)
式(5)中，δRi为信息变化矩阵，δSi为奇异值变化矩阵，Vi为标准正交矩阵，Rmin为子高斯分布的平方根信息矩阵的下界；
Figure GDA0002884026650000041
式(6)中，Qi是标准正交矩阵，Ri是子高斯分布的平方根信息矩阵，δRi为信息变化矩阵，R0为平方根信息矩阵。
3.2.5)根据式(7)计算子高斯分布协方差相比初始高斯分布协方差的减小量；
δPi＝P0-Pi＝R0 -1R0 -T-Ri -1Ri -T＝δYi(δYi)T (7)
式(7)中，δPi为子高斯分布协方差相比初始高斯分布协方差的减小量，δYi为协方差减小量的平方根，δPi半正定且δYi是其矩阵平方根，δYi的函数表达式如式(8)所示；
Figure GDA0002884026650000042
式(8)中，δYi为协方差减小量的平方根，R0为平方根信息矩阵，δRi为信息变化矩阵，矩阵Rdi通过如式(9)所示QR分解得到；
Figure GDA0002884026650000043
式(9)中，Qdi是标准正交矩阵，Rdi是上三角方阵，R0为平方根信息矩阵，δRi为信息变化矩阵，I为n维单位矩阵；
3.2.6)根据初始偏差协方差矩阵P0、子高斯分布协方差相比初始高斯分布协方差的减小量δPi确定子高斯分布协方差矩阵Pi。
优选地，步骤3.3)的详细步骤包括：
3.3.1)记nδY为协方差减小量的平方根矩阵δYi的列数、奇异值变化矩阵δSi的行数，设定打靶次数M作为子高斯分布个数，初始化i；
3.3.2)记当前打靶次数为第i次；
3.3.3)产生nδY维列向量ηi，其各分量服从标准高斯分布，且相互独立；
3.3.4)根据μi＝μ0+δYiηi计算第i个子高斯分布的均值μi，其中μ0为初始状态均值，δYi为协方差减小量的平方根，ηi为nδY维列向量；
3.3.5)确定第i个子高斯分布对应的平方根信息矩阵Ri以及协方差矩阵Pi；
3.3.6)确定第i个子高斯分布对应的权重为wi＝1/M，M为打靶次数；
3.3.7)将变量i加1，判断i小于N是否成立，如果小于N则跳转执行步骤3.3.2)；否则，跳转执行下一步；
3.3.8)产生高斯和形式的初始偏差概率密度函数如式(10)所示；
Figure GDA0002884026650000051
式(10)中，p(t,x0)为高斯和形式的初始偏差概率密度函数，wi为第i个子高斯分布对应的权重，M为打靶次数，pg(x0；μi,Pi)为第i个子高斯分布的概率密度函数，且第i个子高斯分布的概率密度函数pg(x0；μi,Pi)的函数表达式如式(11)所示；
Figure GDA0002884026650000052
式(11)中，n为***维数，Pi为第i个子高斯分布对应的协方差矩阵，x0为初始状态，μi为第i个子高斯分布的均值。
优选地，步骤4)的详细步骤包括：
4.1)根据高斯分布的特性，确定子高斯分布的任意奇数阶中心矩均为0；
4.2)对第i个子高斯分布而言，令
Figure GDA0002884026650000053
其中a＝1,2,...,n，xa是状态变量x的第a个元素，
Figure GDA0002884026650000054
是状态均值μi的第a个元素；根据式(12)计算子高斯分布偏差的二阶中心矩；
E{δxaδxb}＝Pi ab,a,b∈[1,2,...,n] (12)
式(12)中，E{·}表示求括号内物理量的期望值，E{δxaδxb}为子高斯分布偏差的二阶中心矩，δxa为状态变量x与状态均值μi第a个元素的偏差，δxb为状态变量x与状态均值μi第b个元素的偏差，Pi ab表示第i个子高斯分布对应的协方差矩阵Pi的第a行第b列元素；
4.3)在已知子高斯分布偏差的二阶中心矩的基础上，令m初始值为2且依次以2作为增量递增，在第m阶中心矩和协方差矩阵Pi的基础上根据式(13)推导得到子高斯分布的第m+2阶中心矩，直至计算到子高斯分布的前2N阶中心矩，其中N为偏差演化阶数；
Figure GDA0002884026650000055
式(13)中，E{·}表示求括号内物理量的期望值，E{δxaδxbδxc...δxd}为子高斯分布偏差的m+2阶中心矩的第(a,b,c,...,d)个元素，a,b,c,...,d为m+2个对状态变量x的元素的索引，(a,b,c,...,d)为对m+2阶中心矩的元素的索引；Pi ab表示第i个子高斯分布的协方差矩阵Pi的第a行第b列元素，E{δxc...δxd}为子高斯分布偏差的m阶中心矩的第(c,...,d)个元素；Pi ac表示第i个子高斯分布的协方差矩阵Pi的第a行第c列元素，E{δxb...δxd}为子高斯分布偏差的m阶中心矩的第(b,...,d)个元素；以此类推，直至Pi ad表示第i个子高斯分布的协方差矩阵Pi的第a行第d列元素，E{δxbδxc...}为子高斯分布偏差的m阶中心矩的第(b,c,...)个元素。
优选地，步骤5)的详细步骤包括：
5.1)遍历选择一个子高斯分布作为当前的第i个子高斯分布；
5.2)对第i个子高斯分布，根据式(14)将初始状态相对标称值
Figure GDA0002884026650000064
的初始状态偏差δx0变换成第i个子高斯分布的均值μi和相应偏差δxi；
δx0＝μi+δxi (14)
式(14)中，δx0为初始状态偏差，μi为第i个子高斯分布的均值，δxi为第i个子高斯分布的偏差；
5.3)将变换成第i个子高斯分布的均值μi和相应偏差δxi的初始状态偏差δx0代入终端状态关于初始状态偏差δx0的高阶泰勒展开多项式中，得到式(15)所示终端状态；
Figure GDA0002884026650000061
式(15)中，xf为终端状态，T为高阶泰勒展开多项式，T的上标N为偏差演化阶数、T的下标xf为终端状态，δxi为第i个子高斯分布的偏差；
5.4)根据式(16)计算终端时刻第i个子高斯分布的均值μi_f作为式(15)所示终端状态的期望值；
Figure GDA0002884026650000062
式(16)中，μi_f为终端时刻第i个子高斯分布的均值，xf为终端状态，T为高阶泰勒展开多项式，T的上标N为偏差演化阶数、T的下标xf为终端状态，δxi为第i个子高斯分布的偏差；
5.5)根据式(17)将终端时刻第i个子高斯分布的协方差矩阵Pi_f表示成关于初始子高斯分布偏差δxi的2N阶泰勒展开式并求出第i个子高斯分布对应的期望值；
Figure GDA0002884026650000063
式(17)中，Pi_f表示终端时刻第i个子高斯分布的协方差矩阵，xf为终端状态，μi_f为终端时刻第i个子高斯分布的均值，T为高阶泰勒展开多项式，T的上标中的N为偏差演化阶数，δxi为第i个子高斯分布的偏差；
5.6)判断所有的子高斯分布是否已经遍历完毕，如果尚未遍历完毕，则跳转执行步骤5.2)；否则，跳转执行下一步；
5.7)在得到所有子高斯分布在终端时刻的均值μi_f和协方差矩阵Pi_f的基础上，根据式(18)求取终端时刻***状态分布的精确概率密度函数并输出；
Figure GDA0002884026650000071
式(18)中，p(t,xf)表示终端时刻***状态分布的精确概率密度函数，t为终端时刻，xf为终端状态，wi为第i个子高斯分布对应的权重，M为打靶次数，pg(xf；μi_f,Pi_f)为第i个子高斯分布的概率密度函数，μi_f为终端时刻第i个子高斯分布的均值，Pi_f表示终端时刻第i个子高斯分布的协方差矩阵。
本发明基于微分代数与高斯和的非线性***状态偏差演化方法具有下述优点：
1、本发明基于微分代数方法，能够以任意指定阶精度预报非线性***的状态，而且采用了一种递推算法计算任意阶的高斯分布中心矩，因此基于微分代数与高斯和的非线性***状态偏差演化方法可以自动拓展至任意指定阶偏差演化精度，无需手动推导动力学方程的高阶偏导数，具有使用方便、计算精度高的优点，适用于非线性强和长时间预报的非线性***偏差演化分析问题。
2、本发明采用高斯和模型表征***状态偏差分布的概率密度函数，具有对非线性***偏差非高斯性的精确描述能力，而且各子高斯分布预报均基于同一终端状态关于初始状态偏差的高阶泰勒展开多项式，不需为每个子高斯分布分别预报终端状态，因此基于微分代数与高斯和的非线性***状态偏差演化方法在精确描述偏差分布及其非高斯性的同时，仍能保持相对Monte Carlo仿真方法的显著效率优势。
3、在航天器轨道动力学领域，本发明可应用于航天器绝对或相对轨迹预报及其偏差演化分析、轨道确定误差分析等，获得的偏差演化信息可进一步用于空间碰撞预警、航天器轨迹安全性分析与设计等，具有使用方便、计算精度和计算效率高等优点，适用于强非线性***的长时间偏差演化分析问题。
附图说明
图1为本发明实施例一方法的基本流程示意图。
图2为本发明实施例一中求子高斯分布协方差矩阵Pi的流程示意图。
图3为本发明实施例一的终端状态偏差分布对比图。
图4为本发明实施例一的终端概率密度函数对比图。
图5为本发明实施例二的终端状态偏差分布对比图。
图6为本发明实施例二的终端概率密度函数对比图。
具体实施方式
实施例一：
如图1所示，本实施例基于微分代数与高斯和的非线性***状态偏差演化方法步骤包括：
1)输入非线性***的动力学方程、偏差演化阶数N、初始状态相对标称值
Figure GDA0002884026650000082
和初始偏差协方差矩阵P0；
2)采用微分代数方法，针对非线性***的动力学方程预报非线性***的终端状态，并表示为关于初始状态偏差δx0的高阶泰勒展开多项式；
3)针对初始偏差协方差矩阵P0确定子高斯分布协方差矩阵Pi，针对子高斯分布协方差矩阵Pi中的多个子高斯分布进行打靶拟合初始偏差分布生成高斯和模型；
4)针对生成高斯和模型计算给定阶数的子高斯分布高阶中心矩；
5)基于终端状态关于初始状态偏差的高阶多项式和子高斯分布的各阶中心矩，确定各子高斯分布在终端时刻的均值和协方差矩阵，进而给出高斯和形式的终端状态偏差分布概率密度函数。
本实施例基于微分代数与高斯和的非线性***状态偏差演化方法采用微分代数方法，预报非线性***终端状态，并将其表示为关于初始状态偏差的高阶多项式；采用高斯和模型，以多个子高斯分布加权和的形式拟合初始偏差分布；利用子高斯分布的协方差矩阵，计算给定阶数的子高斯分布高阶中心矩；基于终端状态关于初始状态偏差的高阶多项式和子高斯分布的各阶矩，确定各子高斯分布在终端时刻的均值和协方差矩阵。本发明面向一般的非线性动力学***，动力学模型可以是状态转移方程形式，或常微分方程形式，可自动进行任意高阶的状态预报及其偏差演化分析。在航天器轨道动力学领域，本发明可具体应用于航天器绝对或相对轨迹预报及其偏差演化分析、轨道确定误差分析等，获得的偏差演化信息可进一步用于空间碰撞预警、航天器轨迹安全性分析与设计等，具有使用方便、计算精度和计算效率高等优点，适用于强非线性***的长时间偏差演化分析问题。
本实施例基于微分代数与高斯和的非线性***状态偏差演化方法面向一般的非线性动力学***，动力学模型可以是状态转移方程形式，或常微分方程形式，可自动进行任意高阶的状态预报及其偏差演化分析。本实施例中，步骤1)中非线性***的动力学方程为非线性映射形式xf＝f(x0)或常微分方程形式
Figure GDA0002884026650000081
且所述非线性映射形式的函数f(·)或者常微分方程形式的函数g(·)均为任意从n维实数域Rn到Rn的非线性映射。
本实施例中非线性***为近地航天器开普勒轨道，***维数n＝6，状态变量x为地球惯性坐标系下航天器的位置和速度矢量，即x＝[rT vT]T，其中r是位置矢量，v是速度矢量。动力学方程为开普勒轨道解析解，其具体求解形式在现有航天器轨道动力学专著中均可找到，可将航天器终端状态xf看作初始状态x0和轨道预报时间tf的非线性函数，即类似xf＝f(tf,x0)的形式，给定初始时刻航天器经典轨道根数标称值为：
E0＝[a0,e0,i0,Ω0,ω0,f0]＝[6778138,0.001,0.001,0.001,0.001,0.001]
其中a0为半长轴，e0为偏心率，i0为轨道倾角，Ω0为升交点赤经，ω0为近拱点角距，f0为真近点角，长度单位为m，角度单位为rad。根据初始标称轨道根数E0可以计算出初始时刻航天器标称状态
Figure GDA0002884026650000091
给定轨道预报时间tf为一个轨道周期，即tf＝5553.63s。
此外，本实施例给定偏差演化阶数N＝4，给定的初始偏差协方差矩阵P0为：
Figure GDA0002884026650000092
本实施例中，步骤2)的详细步骤包括：
2.1)将初始状态x0初始化为式(1)所示包含初始状态相对标称值
Figure GDA0002884026650000093
的偏差形式；
Figure GDA0002884026650000094
式(1)中，x0为初始状态，
Figure GDA0002884026650000095
为初始状态相对标称值，δx0为初始状态偏差；
2.2)将初始状态x0代入非线性***的动力学方程，应用微分代数方法，求出关于初始状态偏差δx0的高阶泰勒展开多项式描述的终端状态如式(2)所示；
Figure GDA0002884026650000096
式(2)中，xf为终端状态，T为高阶泰勒展开多项式，T的上标N为偏差演化阶数、T的下标xf为高阶泰勒展开多项式描述的内容为终端状态，δx0为初始状态偏差。
本实施例中，步骤3)的详细步骤包括：
3.1)针对初始偏差协方差矩阵P0的逆进行Cholesky分解求平方根信息矩阵R0，且满足P0＝R0 -1R0 -T；
3.2)给定子高斯分布的平方根信息矩阵的下界Rmin，求子高斯分布协方差矩阵Pi；
3.3)将多个子高斯分布打靶生成高斯和模型来拟合初始高斯分布。
如图2所示，步骤3.2)中求子高斯分布协方差矩阵Pi的详细步骤包括：
3.2.1)采用式(3)所示函数表达式计算矩阵
Figure GDA0002884026650000101
的奇异值分解，其中R0为平方根信息矩阵，Rmin为子高斯分布的平方根信息矩阵的下界；
Figure GDA0002884026650000102
式(3)中，Ui和Vi为标准正交矩阵，对角阵Si＝diag(σi1,...,σin)，其中对角线元素σi1,...,σin是正奇异值，R0为平方根信息矩阵，Rmin为子高斯分布的平方根信息矩阵的下界；本实施例中，将子高斯分布平方根信息矩阵的下界设定为Rmin＝5R0；
3.2.2)对任意k＝1,...,n，判断对角阵Si中的第k个对角线元素σik大于或等于1是否成立，如果成立则所求子高斯分布的平方根信息矩阵Ri为平方根信息矩阵R0，跳转执行步骤3.2.5)；否则，跳转执行下一步；
3.2.3)建立式(4)所示n维对角矩阵，删除n维对角矩阵的全0行，得到奇异值变化矩阵δSi；
Figure GDA0002884026650000103
式(4)中，δSi_full为n维对角矩阵，σi1为对角阵Si中的第1个对角线元素，σin为对角阵Si中的第n个对角线元素；
3.2.4)根据式(5)求取信息变化矩阵δRi，并根据式(6)求解子高斯分布的平方根信息矩阵Ri；
δRi＝δSiVi TRmin (5)
式(5)中，δRi为信息变化矩阵，δSi为奇异值变化矩阵，Vi为标准正交矩阵，Rmin为子高斯分布的平方根信息矩阵的下界；
Figure GDA0002884026650000104
式(6)中，Qi是标准正交矩阵，Ri是子高斯分布的平方根信息矩阵，δRi为信息变化矩阵，R0为平方根信息矩阵。
3.2.5)根据式(7)计算子高斯分布协方差相比初始高斯分布协方差的减小量；
δPi＝P0-Pi＝R0 -1R0 -T-Ri -1Ri -T＝δYi(δYi)T (7)
式(7)中，δPi为子高斯分布协方差相比初始高斯分布协方差的减小量，δYi为协方差减小量的平方根，δPi半正定且δYi是其矩阵平方根，δYi的函数表达式如式(8)所示；
Figure GDA0002884026650000111
式(8)中，δYi为协方差减小量的平方根，R0为平方根信息矩阵，δRi为信息变化矩阵，矩阵Rdi通过如式(9)所示QR分解得到；
Figure GDA0002884026650000112
式(9)中，Qdi是标准正交矩阵，Rdi是上三角方阵，R0为平方根信息矩阵，δRi为信息变化矩阵，I为n维单位矩阵；
3.2.6)根据初始偏差协方差矩阵P0、子高斯分布协方差相比初始高斯分布协方差的减小量δPi确定子高斯分布协方差矩阵Pi。
本实施例中，步骤3.3)的详细步骤包括：
3.3.1)记nδY为协方差减小量的平方根矩阵δYi的列数、奇异值变化矩阵δSi的行数，设定打靶次数M作为子高斯分布个数，初始化i；本实施例中，子高斯分布数量M＝5000；
3.3.2)记当前打靶次数为第i次；
3.3.3)产生nδY维列向量ηi，其各分量服从标准高斯分布，且相互独立；
3.3.4)根据μi＝μ0+δYiηi计算第i个子高斯分布的均值μi，其中μ0为初始状态均值，δYi为协方差减小量的平方根，ηi为nδY维列向量；
3.3.5)确定第i个子高斯分布对应的平方根信息矩阵Ri以及协方差矩阵Pi；
3.3.6)确定第i个子高斯分布对应的权重为wi＝1/M，M为打靶次数；
3.3.7)将变量i加1，判断i小于N是否成立，如果小于N则跳转执行步骤3.3.2)；否则，跳转执行下一步；
3.3.8)产生高斯和形式的初始偏差概率密度函数如式(10)所示；
Figure GDA0002884026650000113
式(10)中，p(t,x0)为高斯和形式的初始偏差概率密度函数，wi为第i个子高斯分布对应的权重，M为打靶次数，pg(x0；μi,Pi)为第i个子高斯分布的概率密度函数，且第i个子高斯分布的概率密度函数pg(x0；μi,Pi)的函数表达式如式(11)所示；
Figure GDA0002884026650000121
式(11)中，n为***维数，Pi为第i个子高斯分布对应的协方差矩阵，x0为初始状态，μi为第i个子高斯分布的均值。
本实施例中，步骤4)的详细步骤包括：
4.1)根据高斯分布的特性，确定子高斯分布的任意奇数阶中心矩均为0；
例如，子高斯分布的一阶中心矩为：
E{δxa}＝0
其中，E{·}表示求括号内物理量的期望值，E{δxa}为子高斯分布偏差的一阶中心矩，xa是状态变量x的第a个元素，δxa为状态变量x与状态均值μi第a个元素的偏差。
例如，子高斯分布的三阶中心矩为：
E{δxaδxbδxc}＝0,a,b,c∈[1,2,...,n]
其中，E{·}表示求括号内物理量的期望值，E{δxaδxbδxc}为子高斯分布偏差的一阶中心矩，xa是状态变量x的第a个元素，δxa为状态变量x与状态均值μi第a个元素的偏差，xb为状态变量x的第b个元素，δxb为状态变量x与状态均值μi第b个元素的偏差，xc为状态变量x的第c个元素，δxc为状态变量x与状态均值μi第c个元素的偏差。更高奇数阶中心矩类似，均为0。
4.2)对第i个子高斯分布而言，令
Figure GDA0002884026650000122
其中a＝1,2,...,n，xa是状态变量x的第a个元素，
Figure GDA0002884026650000123
是状态均值μi的第a个元素；根据式(12)计算子高斯分布偏差的二阶中心矩；
E{δxaδxb}＝Pi ab,a,b∈[1,2,...,n] (12)
式(12)中，E{·}表示求括号内物理量的期望值，E{δxaδxb}为子高斯分布偏差的二阶中心矩，δxa为状态变量x与状态均值μi第a个元素的偏差，δxb为状态变量x与状态均值μi第b个元素的偏差，Pi ab表示第i个子高斯分布对应的协方差矩阵Pi的第a行第b列元素；
4.3)在已知子高斯分布偏差的二阶中心矩的基础上，令m初始值为2且依次以2作为增量递增(m＝2,4,6,...)，在第m阶中心矩和协方差矩阵Pi的基础上根据式(13)推导得到子高斯分布的第m+2阶中心矩，直至计算到子高斯分布的前2N阶中心矩，其中N为偏差演化阶数；
Figure GDA0002884026650000131
式(13)中，E{·}表示求括号内物理量的期望值，E{δxaδxbδxc...δxd}为子高斯分布偏差的m+2阶中心矩的第(a,b,c,...,d)个元素，a,b,c,...,d为m+2个对状态变量x的元素的索引，(a,b,c,...,d)为对m+2阶中心矩的元素的索引；Pi ab表示第i个子高斯分布的协方差矩阵Pi的第a行第b列元素，E{δxc...δxd}为子高斯分布偏差的m阶中心矩的第(c,...,d)个元素；Pi ac表示第i个子高斯分布的协方差矩阵Pi的第a行第c列元素，E{δxb...δxd}为子高斯分布偏差的m阶中心矩的第(b,...,d)个元素；以此类推，直至Pi ad表示第i个子高斯分布的协方差矩阵Pi的第a行第d列元素，E{δxbδxc...}为子高斯分布偏差的m阶中心矩的第(b,c,...)个元素。
本实施例中，步骤5)的详细步骤包括：
5.1)遍历选择一个子高斯分布作为当前的第i个子高斯分布；
5.2)对第i个子高斯分布，根据式(14)将初始状态相对标称值
Figure GDA0002884026650000133
的初始状态偏差δx0变换成第i个子高斯分布的均值μi和相应偏差δxi；
δx0＝μi+δxi (14)
式(14)中，δx0为初始状态偏差，μi为第i个子高斯分布的均值，δxi为第i个子高斯分布的偏差；
5.3)将变换成第i个子高斯分布的均值μi和相应偏差δxi的初始状态偏差δx0代入终端状态关于初始状态偏差δx0的高阶泰勒展开多项式中，得到式(15)所示终端状态；
Figure GDA0002884026650000132
式(15)中，xf为终端状态，T为高阶泰勒展开多项式，T的上标N为偏差演化阶数、T的下标xf为终端状态，δxi为第i个子高斯分布的偏差；
5.4)根据式(16)计算终端时刻第i个子高斯分布的均值μi_f作为式(15)所示终端状态的期望值；
Figure GDA0002884026650000141
式(16)中，μi_f为终端时刻第i个子高斯分布的均值，xf为终端状态，T为高阶泰勒展开多项式，T的上标N为偏差演化阶数、T的下标xf为终端状态，δxi为第i个子高斯分布的偏差；
5.5)根据式(17)将终端时刻第i个子高斯分布的协方差矩阵Pi_f表示成关于初始子高斯分布偏差δxi的2N阶泰勒展开式并求出第i个子高斯分布对应的期望值；
Figure GDA0002884026650000142
式(17)中，Pi_f表示终端时刻第i个子高斯分布的协方差矩阵，xf为终端状态，μi_f为终端时刻第i个子高斯分布的均值，T为高阶泰勒展开多项式，T的上标中的N为偏差演化阶数，δxi为第i个子高斯分布的偏差；
5.6)判断所有的子高斯分布是否已经遍历完毕，如果尚未遍历完毕，则跳转执行步骤5.2)；否则，跳转执行下一步；
5.7)在得到所有子高斯分布在终端时刻的均值μi_f和协方差矩阵Pi_f的基础上，根据式(18)求取终端时刻***状态分布的精确概率密度函数并输出；
Figure GDA0002884026650000143
式(18)中，p(t,xf)表示终端时刻***状态分布的精确概率密度函数，t为终端时刻，xf为终端状态，wi为第i个子高斯分布对应的权重，M为打靶次数，pg(xf；μi_f,Pi_f)为第i个子高斯分布的概率密度函数，μi_f为终端时刻第i个子高斯分布的均值，Pi_f表示终端时刻第i个子高斯分布的协方差矩阵。
参见图3所示终端状态偏差分布对比图中，子高斯分布表示本发明预报得到的终端时刻子高斯偏差分布情况，MonteCarlo打靶点表示传统MonteCarlo仿真方法得到的终端时刻打靶点分布情况，MonteCarlo协方差表示根据MonteCarlo打靶点统计出的协方差3-sigma椭球，MonteCarlo均值表示根据MonteCarlo打靶点统计出的分布均值，线性化协方差表示传统线性化方法得到的终端时刻协方差3-sigma椭球，线性化均值表示传统线性化方法得到的终端时刻分布均值。参见图3可知，与精确的MonteCarlo仿真相比，由传统线性化方法预报的终端均值和协方差矩阵均存在很大误差，而本实施例基于微分代数与高斯和的非线性***状态偏差演化方法用高斯和模型很好地拟合了MonteCarlo仿真样本点的分布，验证了本方法的正确性和有效性。
参见图4所示终端概率密度函数对比图中，高斯和表示本发明得到的终端高斯和模型概率密度函数，MonteCarlo表示根据MonteCarlo打靶点统计出的概率密度函数，线性化方法表示传统线性化方法得到的终端时刻概率密度函数。进一步参见图4可知，本实施例基于微分代数与高斯和的非线性***状态偏差演化方法计算出的概率密度函数与MonteCarlo仿真统计结果一致，而且明显概率密度函数已不再呈高斯分布形式，同时可发现传统线性化方法则存在明显的误差。因此说明本方法可以精确预报强非线性***的非高斯偏差演化情况。
综上所述，本实施例基于微分代数与高斯和的非线性***状态偏差演化方法为一种高精度高效率预报非线性***状态偏差演化过程的方法，基于微分代数预报***终端状态，应用一种递推算法计算高斯分布的任意阶中心矩，可以自动拓展至任意阶精度的偏差演化预报，并且采用高斯和模型描述状态偏差的非高斯分布特性，在精确描述偏差分布及其非高斯性的同时仍能保持相对MonteCarlo仿真方法的显著效率优势，适用于强非线性***和长时间精确预报问题。其中，动力学方程可以是任意非线性映射形式或非线性常微分方程形式，本实施例的动力学方程是非线性映射形式的开普勒轨道解析解。本实施例基于微分代数与高斯和的非线性***状态偏差演化方法面向一般非线性动力学***，基于微分代数和高斯和方法的，可自动拓展至任意高阶精度的，实现通用非线性***状态偏差精确演化。
实施例二：
本实施例的过程与实施例一基本相同，其主要不同点为：
步骤2.1)中，非线性***为地球静止轨道考虑太阳光压的航天器非线性相对运动轨迹，***维数n＝6，状态变量x为从航天器在主航天器轨道坐标系下的相对位置和速度矢量，即x＝[rT vT]T，其中r＝[x,y,z]T是相对位置矢量，v＝[vx,vy,vz]T是相对速度矢量。
本实施例中，非线性***的动力学方程为常微分方程形式，终端状态通过变步长龙格库塔积分方法得到。
本实施例中，非线性***的动力学方程的具体表达式为：
Figure GDA0002884026650000161
其中，a为主航天器轨道半长轴，本实施例中即为地球静止轨道半径a＝42164200m，静止轨道航天器的轨道角速度为n＝7.2921×10-5rad/s，地球引力常数为μ＝3.986×1014m3/s2，太阳光压摄动加速度为γ＝[γx,γy,γz]T，其具体表达形式可参考文献：Fehse,W.,“Disturbance by Solar Pressure of Rendezvous Trajectories inGEO,”5th ICATT,Noordwijk,Netherlands,May 2012。
给定初始时刻航天器经典轨道根数标称值为：
Figure GDA0002884026650000162
给定轨道预报时间tf为一个轨道周期，本实施例中即为tf＝86164s。
步骤2.1)中，给定初始状态协方差矩阵P0为：
Figure GDA0002884026650000163
相比实施例一而言，本实施例非线性***变为地球静止轨道考虑太阳光压的航天器非线性相对运动轨迹，动力学方程变为常微分方程形式，终端状态通过变步长龙格库塔积分方法得到。
参见图5所示终端状态偏差分布对比图中，子高斯分布表示本发明预报得到的终端时刻子高斯偏差分布情况，MonteCarlo打靶点表示传统MonteCarlo仿真方法得到的终端时刻打靶点分布情况，MonteCarlo协方差表示根据MonteCarlo打靶点统计出的协方差3-sigma椭球，MonteCarlo均值表示根据MonteCarlo打靶点统计出的分布均值，线性化协方差表示传统线性化方法得到的终端时刻协方差3-sigma椭球，线性化均值表示传统线性化方法得到的终端时刻分布均值。参见图6所示终端概率密度函数对比图中，高斯和表示本发明得到的终端高斯和模型概率密度函数，MonteCarlo表示根据MonteCarlo打靶点统计出的概率密度函数，线性化方法表示传统线性化方法得到的终端时刻概率密度函数。参见图5所示终端状态偏差分布对比和图6所示终端概率密度函数对比，本实施例基于微分代数与高斯和的非线性***状态偏差演化方法仍然能很好地拟合了MonteCarlo仿真结果，有效捕获航天器相对轨迹偏差分布的非高斯性，并准确预报终端时刻偏差分布的概率密度函数。本实施例验证了本实施例基于微分代数与高斯和的非线性***状态偏差演化方法的普适性，对一般的非线性***，动力学方程可以是非线性映射形式或常微分方程形式，均能精确预报其状态偏差的演化情况。
以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (6)

1.一种基于微分代数与高斯和的非线性***状态偏差演化方法，用于航天器轨迹偏差演化，其特征为该方法的步骤包括：
1)输入非线性***的动力学方程、偏差演化阶数N、初始状态相对标称值
Figure FDA0002884026640000011
和初始偏差协方差矩阵P0；
2)采用微分代数方法，针对非线性***的动力学方程预报非线性***的终端状态，并表示为关于初始状态偏差δx0的高阶泰勒展开多项式；
3)针对初始偏差协方差矩阵P0确定子高斯分布协方差矩阵Pi，针对子高斯分布协方差矩阵Pi中的多个子高斯分布进行打靶拟合初始偏差分布生成高斯和模型；
4)针对生成高斯和模型计算给定阶数的子高斯分布高阶中心矩；
5)基于终端状态关于初始状态偏差的高阶多项式和子高斯分布的各阶中心矩，确定各子高斯分布在终端时刻的均值和协方差矩阵，进而给出高斯和形式的终端状态偏差分布概率密度函数；
步骤4)的详细步骤包括：
4.1)根据高斯分布的特性，确定子高斯分布的任意奇数阶中心矩均为0；
4.2)对第i个子高斯分布而言，令
Figure FDA0002884026640000012
其中a＝1,2,...,n，xa是状态变量x的第a个元素，
Figure FDA0002884026640000013
是状态均值μi的第a个元素；根据式(12)计算子高斯分布偏差的二阶中心矩；
E{δxaδxb}＝Pi ab,a,b∈[1,2,...,n] (12)
式(12)中，E{·}表示求括号内物理量的期望值，E{δxaδxb}为子高斯分布偏差的二阶中心矩，δxa为状态变量x与状态均值μi第a个元素的偏差，δxb为状态变量x与状态均值μi第b个元素的偏差，Pi ab表示第i个子高斯分布对应的协方差矩阵Pi的第a行第b列元素；
4.3)在已知子高斯分布偏差的二阶中心矩的基础上，令m初始值为2且依次以2作为增量递增，在第m阶中心矩和协方差矩阵Pi的基础上根据式(13)推导得到子高斯分布的第m+2阶中心矩，直至计算到子高斯分布的前2N阶中心矩，其中N为偏差演化阶数；
Figure FDA0002884026640000014
式(13)中，E{·}表示求括号内物理量的期望值，E{δxaδxbδxc...δxd}为子高斯分布偏差的m+2阶中心矩的第(a,b,c,...,d)个元素，(a,b,c,...,d)为对m+2阶中心矩的元素的索引；Pi ab表示第i个子高斯分布的协方差矩阵Pi的第a行第b列元素，E{δxc...δxd}为子高斯分布偏差的m阶中心矩的第(c,...,d)个元素；Pi ac表示第i个子高斯分布的协方差矩阵Pi的第a行第c列元素，E{δxb...δxd}为子高斯分布偏差的m阶中心矩的第(b,...,d)个元素；以此类推，直至Pi ad表示第i个子高斯分布的协方差矩阵Pi的第a行第d列元素，E{δxbδxc...}为子高斯分布偏差的m阶中心矩的第(b,c,...)个元素；
步骤5)的详细步骤包括：
5.1)遍历选择一个子高斯分布作为当前的第i个子高斯分布；
5.2)对第i个子高斯分布，根据式(14)将初始状态相对标称值
Figure FDA0002884026640000023
的初始状态偏差δx0变换成第i个子高斯分布的均值μi和相应偏差δxi；
δx0＝μi+δxi (14)
式(14)中，δx0为初始状态偏差，μi为第i个子高斯分布的均值，δxi为第i个子高斯分布的偏差；
5.3)将变换成第i个子高斯分布的均值μi和相应偏差δxi的初始状态偏差δx0代入终端状态关于初始状态偏差δx0的高阶泰勒展开多项式中，得到式(15)所示终端状态；
Figure FDA0002884026640000021
式(15)中，xf为终端状态，T为高阶泰勒展开多项式，T的上标N为偏差演化阶数、T的下标xf为终端状态，δxi为第i个子高斯分布的偏差；
5.4)根据式(16)计算终端时刻第i个子高斯分布的均值μi_f作为式(15)所示终端状态的期望值；
Figure FDA0002884026640000022
式(16)中，μi_f为终端时刻第i个子高斯分布的均值，xf为终端状态，T为高阶泰勒展开多项式，T的上标N为偏差演化阶数、T的下标xf为终端状态，δxi为第i个子高斯分布的偏差；
5.5)根据式(17)将终端时刻第i个子高斯分布的协方差矩阵Pi_f表示成关于初始子高斯分布偏差δxi的2N阶泰勒展开式并求出第i个子高斯分布对应的期望值；
Figure FDA0002884026640000031
式(17)中，Pi_f表示终端时刻第i个子高斯分布的协方差矩阵，xf为终端状态，μi_f为终端时刻第i个子高斯分布的均值，T为高阶泰勒展开多项式，T的上标中的N为偏差演化阶数，δxi为第i个子高斯分布的偏差；
5.6)判断所有的子高斯分布是否已经遍历完毕，如果尚未遍历完毕，则跳转执行步骤5.2)；否则，跳转执行下一步；
5.7)在得到所有子高斯分布在终端时刻的均值μi_f和协方差矩阵Pi_f的基础上，根据式(18)求取终端时刻***状态分布的精确概率密度函数并输出；
Figure FDA0002884026640000032
式(18)中，p(t,xf)表示终端时刻***状态分布的精确概率密度函数，t为终端时刻，xf为终端状态，wi为第i个子高斯分布对应的权重，M为打靶次数，pg(xf；μi_f,Pi_f)为第i个子高斯分布的概率密度函数，μi_f为终端时刻第i个子高斯分布的均值，Pi_f表示终端时刻第i个子高斯分布的协方差矩阵。
2.根据权利要求1所述的基于微分代数与高斯和的非线性***状态偏差演化方法，其特征为，步骤1)中非线性***的动力学方程为非线性映射形式xf＝f(x0)或常微分方程形式
Figure FDA0002884026640000033
且所述非线性映射形式的函数f(·)或者常微分方程形式的函数g(·)均为任意从n维实数域Rn到Rn的非线性映射，其中x0为初始状态，xf为终端状态，x为状态变量。
3.根据权利要求1所述的基于微分代数与高斯和的非线性***状态偏差演化方法，其特征为，步骤2)的详细步骤包括：
2.1)将初始状态x0初始化为式(1)所示包含初始状态相对标称值
Figure FDA0002884026640000034
的偏差形式；
Figure FDA0002884026640000035
式(1)中，x0为初始状态，
Figure FDA0002884026640000036
为初始状态相对标称值，δx0为初始状态偏差；
2.2)将初始状态x0代入非线性***的动力学方程，应用微分代数方法，求出关于初始状态偏差δx0的高阶泰勒展开多项式描述的终端状态如式(2)所示；
Figure FDA0002884026640000037
式(2)中，xf为终端状态，T为高阶泰勒展开多项式，T的上标N为偏差演化阶数、T的下标xf为高阶泰勒展开多项式描述的内容为终端状态，δx0为初始状态偏差。
4.根据权利要求1所述的基于微分代数与高斯和的非线性***状态偏差演化方法，其特征为，步骤3)的详细步骤包括：
3.1)针对初始偏差协方差矩阵P0的逆进行Cholesky分解求平方根信息矩阵R0，且满足P0＝R0 -1R0 -T；
3.2)给定子高斯分布的平方根信息矩阵的下界Rmin，求子高斯分布协方差矩阵Pi；
3.3)将多个子高斯分布打靶生成高斯和模型来拟合初始高斯分布。
5.根据权利要求4所述的基于微分代数与高斯和的非线性***状态偏差演化方法，其特征为，步骤3.2)中求子高斯分布协方差矩阵Pi的详细步骤包括：
3.2.1)采用式(3)所示函数表达式计算矩阵
Figure FDA0002884026640000041
的奇异值分解，其中R0为平方根信息矩阵，Rmin为子高斯分布的平方根信息矩阵的下界；
Figure FDA0002884026640000042
式(3)中，Ui和Vi为标准正交矩阵，对角阵Si＝diag(σi1,...,σin)，其中对角线元素σi1,...,σin是正奇异值，R0为平方根信息矩阵，Rmin为子高斯分布的平方根信息矩阵的下界；
3.2.2)对任意k＝1,...,n，判断对角阵Si中的第k个对角线元素σik大于或等于1是否成立，如果成立则所求子高斯分布的平方根信息矩阵Ri为平方根信息矩阵R0，跳转执行步骤3.2.5)；否则，跳转执行下一步；
3.2.3)建立式(4)所示n维对角矩阵，删除n维对角矩阵的全0行，得到奇异值变化矩阵δSi；
Figure FDA0002884026640000043
式(4)中，δSi_full为n维对角矩阵，σi1为对角阵Si中的第1个对角线元素，σin为对角阵Si中的第n个对角线元素；
3.2.4)根据式(5)求取信息变化矩阵δRi，并根据式(6)求解子高斯分布的平方根信息矩阵Ri；
δRi＝δSiVi TRmin (5)
式(5)中，δRi为信息变化矩阵，δSi为奇异值变化矩阵，Vi为标准正交矩阵，Rmin为子高斯分布的平方根信息矩阵的下界；
Figure FDA0002884026640000051
式(6)中，Qi是标准正交矩阵，Ri是子高斯分布的平方根信息矩阵，δRi为信息变化矩阵，R0为平方根信息矩阵；
3.2.5)根据式(7)计算子高斯分布协方差相比初始高斯分布协方差的减小量；
δPi＝P0-Pi＝R0 -1R0 -T-Ri -1Ri -T＝δYi(δYi)T (7)
式(7)中，δPi为子高斯分布协方差相比初始高斯分布协方差的减小量，δYi为协方差减小量的平方根，δPi半正定且δYi是δPi的矩阵平方根，δYi的函数表达式如式(8)所示；
Figure FDA0002884026640000052
式(8)中，δYi为协方差减小量的平方根，R0为平方根信息矩阵，δRi为信息变化矩阵，矩阵Rdi通过如式(9)所示QR分解得到；
Figure FDA0002884026640000053
式(9)中，Qdi是标准正交矩阵，Rdi是上三角方阵，R0为平方根信息矩阵，δRi为信息变化矩阵，I为n维单位矩阵；
3.2.6)根据初始偏差协方差矩阵P0、子高斯分布协方差相比初始高斯分布协方差的减小量δPi确定子高斯分布协方差矩阵Pi。
6.根据权利要求5所述的基于微分代数与高斯和的非线性***状态偏差演化方法，其特征为，步骤3.3)的详细步骤包括：
3.3.1)记nδY为协方差减小量的平方根矩阵δYi的列数、奇异值变化矩阵δSi的行数，设定打靶次数M作为子高斯分布个数，初始化i；
3.3.2)记当前打靶次数为第i次；
3.3.3)产生nδY维列向量ηi，其各分量服从标准高斯分布，且相互独立；
3.3.4)根据μi＝μ0+δYiηi计算第i个子高斯分布的均值μi，其中μ0为初始状态均值，δYi为协方差减小量的平方根，ηi为nδY维列向量；
3.3.5)确定第i个子高斯分布对应的平方根信息矩阵Ri以及协方差矩阵Pi；
3.3.6)确定第i个子高斯分布对应的权重为wi＝1/M，M为打靶次数；
3.3.7)将变量i加1，判断i小于N是否成立，如果小于N则跳转执行步骤3.3.2)；否则，跳转执行下一步；
3.3.8)产生高斯和形式的初始偏差概率密度函数如式(10)所示；
Figure FDA0002884026640000061
式(10)中，p(t,x0)为高斯和形式的初始偏差概率密度函数，wi为第i个子高斯分布对应的权重，M为打靶次数，pg(x0；μi,Pi)为第i个子高斯分布的概率密度函数，且第i个子高斯分布的概率密度函数pg(x0；μi,Pi)的函数表达式如式(11)所示；
Figure FDA0002884026640000062
式(11)中，n为***维数，Pi为第i个子高斯分布对应的协方差矩阵，x0为初始状态，μi为第i个子高斯分布的均值。
CN201710551107.9A 2017-07-07 2017-07-07 基于微分代数与高斯和的非线性***状态偏差演化方法 Active CN107402903B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710551107.9A CN107402903B (zh) 2017-07-07 2017-07-07 基于微分代数与高斯和的非线性***状态偏差演化方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710551107.9A CN107402903B (zh) 2017-07-07 2017-07-07 基于微分代数与高斯和的非线性***状态偏差演化方法

Publications (2)

Publication Number Publication Date
CN107402903A CN107402903A (zh) 2017-11-28
CN107402903B true CN107402903B (zh) 2021-02-26

Family

ID=60404964

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710551107.9A Active CN107402903B (zh) 2017-07-07 2017-07-07 基于微分代数与高斯和的非线性***状态偏差演化方法

Country Status (1)

Country Link
CN (1) CN107402903B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108490966B (zh) * 2018-01-31 2021-02-05 中国人民解放军国防科技大学 基于微分代数的静止轨道摄动相对轨迹高阶制导方法
CN109032176B (zh) * 2018-07-25 2021-06-22 西北工业大学 一种基于微分代数的地球同步轨道确定和参数确定方法
CN109508482A (zh) * 2018-10-26 2019-03-22 天津大学 一种用于复杂曲面表面轮廓度误差不确定度的计算方法
CN110442831B (zh) * 2019-07-31 2023-03-24 中国人民解放军国防科技大学 基于非线性偏差演化的空间非合作目标天基搜索方法
CN111899297B (zh) * 2020-08-06 2024-01-23 中国铁建重工集团股份有限公司 线结构光条纹中心提取方法
CN113836697B (zh) * 2021-08-26 2024-01-23 北京理工大学 基于动态混合高斯的空间引力波探测构型稳定性演化方法
CN114329887B (zh) * 2021-11-22 2024-07-12 北京理工大学 一种数据驱动的动态混合高斯航天器轨道误差演化方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106697333A (zh) * 2017-01-12 2017-05-24 北京理工大学 一种航天器轨道控制策略的鲁棒性分析方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101770024B (zh) * 2010-01-25 2011-08-31 上海交通大学 多目标跟踪方法
CN102040008A (zh) * 2010-12-13 2011-05-04 北京航空航天大学 一种用于编队卫星在轨运行安全的防碰控制方法
CN105303052A (zh) * 2015-11-11 2016-02-03 中国人民解放军国防科学技术大学 一种低速接近航天器轨迹安全评价方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106697333A (zh) * 2017-01-12 2017-05-24 北京理工大学 一种航天器轨道控制策略的鲁棒性分析方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Nonlinear Mapping of Uncertainties in Celestial Mechanics;M. Valli et al.;《Journal of Guidance Control and Dynamics》;20130228;第36卷(第1期);第48-63页 *
航天器轨迹偏差的非线性非高斯预报方法;孙振江 等;《宇航学报》;20161130;第37卷(第11期);第1289-1297页 *

Also Published As

Publication number Publication date
CN107402903A (zh) 2017-11-28

Similar Documents

Publication Publication Date Title
CN107402903B (zh) 基于微分代数与高斯和的非线性***状态偏差演化方法
da Silva et al. Ensemble-based state estimator for aerodynamic flows
CN106970643B (zh) 一种解析的卫星非线性相对运动偏差传播分析方法
Ning et al. An implicit UKF for satellite stellar refraction navigation system
CN110378012B (zh) 一种考虑高阶重力场的严格回归轨道设计方法、***及介质
Dyniewicz Space–time finite element approach to general description of a moving inertial load
CN102878997A (zh) 一种大偏心率轨道的星上快速高精度外推方法
Stern et al. Estimation of dynamic stability coefficients for aerodynamic decelerators using CFD
CN114936471A (zh) 一种基于并行计算的航天器碰撞预警分层快速筛选方法
Thompson et al. Aerodynamic moment model calibration from distributed pressure arrays
Huang et al. Reliability-based trajectory optimization using nonintrusive polynomial chaos for Mars entry mission
Schweizer et al. Solving differential-algebraic equation systems: alternative index-2 and index-1 approaches for constrained mechanical systems
CN109709805B (zh) 一种考虑不确定性因素的航天器鲁棒交会轨迹设计方法
Xue et al. A comparison of several nonlinear filters for mobile robot pose estimation
Brock et al. Dynamic CFD simulations of the supersonic inflatable aerodynamic decelerator (SAID) ballistic range tests
CN110231619B (zh) 基于恩克法的雷达交接时刻预报方法及装置
Murman et al. Dynamic analysis of atmospheric-entry probes and capsules
Morris et al. Real-time system identification of quadrotor dynamics
CN111428912A (zh) 一种基于支持向量机的火星探测器轨道预测方法及***
CN114462256B (zh) 非合作小推力机动目标轨道确定方法、装置、设备和介质
Aeschliman et al. A proposed methodology for computational fluid dynamics code verification, calibration, and validation
CN111259337A (zh) 一种基于统计的重残骸实时落点预报方法
Wang et al. Dynamic models of satellite relative motion around an oblate earth
Zhou et al. Inverse simulation system for manual-controlled rendezvous and docking based on artificial neural network
Huang et al. Uncertainty analysis of reachable set for planetary entry using polynomial chaos

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