基于微分代数与高斯和的非线性系统状态偏差演化方法与流程

本发明涉及航天器轨道动力学领域，具体涉及一种基于微分代数与高斯和的非线性系统状态偏差演化方法，用于非线性系统的状态及其偏差预报。

背景技术：

航天器轨迹偏差演化技术是开展空间碰撞预警、航天器轨道确定、导航滤波以及航天器轨迹安全性分析与设计等航天任务的重要基础技术。轨迹偏差的预报精度直接影响了预警、定轨和导航结果的可信度，甚至决定着航天任务的成败。由于航天器轨道动力学模型的实质就是一个常微分方程组或状态转移方程形式的非线性系统，因此航天器轨迹偏差演化问题实质就是非线性系统状态偏差预报问题。

经典的非线性系统偏差演化方法是线性化协方差分析方法和montecarlo仿真方法。线性化协方差分析法，先将非线性动力学系统线性化，再采用线性系统协方差分析方法进行偏差预报计算，该方法使用简单、计算量小，但是对强非线性系统或长时间偏差演化预报问题计算误差较大。montecarlo仿真方法，通过统计多次打靶试验结果可以得到高精度的状态偏差演化结果，但是使用该方法需要进行大量抽样仿真，计算量大，很多情况下计算效率难以满足任务需求。

微分代数方法是一种在数值环境下计算非线性映射关于自变量的高阶偏导数的方法，利用这些高阶偏导数可以对非线性系统进行任意阶次的泰勒展开，进而使用高阶多项式预报系统状态变量。该方法最初由美国密歇根大学粒子物理学教授martinberz提出，基本原理可见参考文献：berzm.differentialalgebraicdescriptionofbeamdynamicstoveryhighorders[j].part.accel.,1988,24(ssc-152):109-124。目前国内尚未有应用该方法进行非线性动力学系统偏差预报的研究或应用成果。

高斯和模型是以多个子高斯分布加权和的形式描述系统状态偏差分布的方法，可用于准确描述非高斯分布的状态偏差概率密度函数。本发明面向一般的非线性动力学系统，系统动力学模型可以是状态转移方程形式，或常微分方程形式，利用微分代数方法求解终端状态关于初始状态偏差的高阶多项式，进而采用高斯和模型拟合初始状态偏差分布，并基于前述高阶多项式预报终端时刻的各子高斯分布，以此来准确描述终端状态的非高斯分布特性。

技术实现要素：

本发明要解决的技术问题：针对现有技术的上述问题，提供一种能够自动拓展至任意指定阶偏差演化精度，无需手动推导动力学方程的高阶偏导数，适用于非线性强和长时间预报的非线性系统偏差演化分析问题，在精确描述偏差分布及其非高斯性的同时，仍能保持相对montecarlo仿真方法的显著效率优势，使用方便、计算精度高的基于微分代数与高斯和的非线性系统状态偏差演化方法。

为了解决上述技术问题，本发明采用的技术方案为：

一种基于微分代数与高斯和的非线性系统状态偏差演化方法，步骤包括：

1)输入非线性系统的动力学方程、偏差演化阶数n、初始状态相对标称值和初始偏差协方差矩阵p0；

2)采用微分代数方法，针对非线性系统的动力学方程预报非线性系统的终端状态，并表示为关于初始状态偏差δx0的高阶泰勒展开多项式；

3)针对初始偏差协方差矩阵p0确定子高斯分布协方差矩阵pi，针对子高斯分布协方差矩阵pi中的多个子高斯分布进行打靶拟合初始偏差分布生成高斯和模型；

4)针对生成高斯和模型计算给定阶数的子高斯分布高阶中心矩；

5)基于终端状态关于初始状态偏差的高阶多项式和子高斯分布的各阶中心矩，确定各子高斯分布在终端时刻的均值和协方差矩阵，进而给出高斯和形式的终端状态偏差分布概率密度函数。

优选地，步骤1)中非线性系统的动力学方程为非线性映射形式xf＝f(x0)或常微分方程形式且所述非线性映射形式的函数f(·)或者常微分方程形式的函数g(·)均为任意从n维实数域rⁿ到rⁿ的非线性映射。

优选地，步骤2)的详细步骤包括：

2.1)将初始状态x0初始化为式(1)所示包含初始状态相对标称值的偏差形式；

式(1)中，x0为初始状态，为初始状态相对标称值，δx0为初始偏差；

2.2)将初始状态x0代入非线性系统的动力学方程，应用微分代数方法，求而出关于初始偏差δx0的高阶泰勒展开多项式描述的终端状态如式(2)所示；

式(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)所示函数表达式计算矩阵的奇异值分解，其中r0为初始偏差协方差矩阵，rmin为子高斯分布的平方根信息矩阵的下界；

式(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；

式(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为子高斯分布的平方根信息矩阵的下界；

式(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)所示；

式(8)中，δyi为协方差减小量的平方根，r0为初始偏差协方差矩阵，δri为信息变化矩阵，矩阵rdi通过如式(9)所示qr分解得到；

式(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)所示；

式(10)中，p(t,x0)为高斯和形式的初始偏差概率密度函数，wi为第i个子高斯分布对应的权重，m为打靶次数，pg(x0；μi,pi)为第i个子高斯分布的概率密度函数，且第i个子高斯分布的概率密度函数pg(x0；μi,pi)的函数表达式如式(11)所示；

式(11)中，n为系统维数，pi为第i个子高斯分布对应的协方差矩阵，x0为初始状态，μi为第i个子高斯分布的均值。

优选地，步骤4)的详细步骤包括：

4.1)根据高斯分布的特性，确定子高斯分布的任意奇数阶中心矩均为0；

4.2)对第i个子高斯分布而言，令其中a＝1,2,...,n，x^a是状态变量x的第a个元素，是状态均值μi的第a个元素；根据式(12)计算子高斯分布偏差的二阶中心矩；

e{δx^aδx^b}＝pi^ab,a,b∈[1,2,...,n](12)

式(12)中，e{·}表示求括号内物理量的期望值，e{δx^aδx^b}为子高斯分布偏差的二阶中心矩，δx^a为状态变量x与状态均值μi第a个元素的偏差，δx^b为状态变量x与状态均值μi第b个元素的偏差，pi^ab表示第i个子高斯分布对应的协方差矩阵pi的第a行第b列元素；

4.3)在已知子高斯分布偏差的二阶中心矩的基础上，令m初始值为2且依次以2作为增量递增，在第m阶中心矩和协方差矩阵pi的基础上根据式(13)推导得到子高斯分布的第m+2阶中心矩，直至计算到子高斯分布的前2n阶中心矩，其中n为偏差演化阶数；

式(13)中，e{·}表示求括号内物理量的期望值，e{δx^aδx^bδx^c...δx^d}为子高斯分布偏差的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{δx^c...δx^d}为子高斯分布偏差的m阶中心矩的第(c,…,d)个元素；pi^ac表示第i个子高斯分布的协方差矩阵pi的第a行第c列元素，e{δx^b...δx^d}为子高斯分布偏差的m阶中心矩的第(b,…,d)个元素；以此类推，直至pi^ad表示第i个子高斯分布的协方差矩阵pi的第a行第d列元素，e{δx^bδx^c...}为子高斯分布偏差的m阶中心矩的第(b,c,…)个元素。

优选地，步骤5)的详细步骤包括：

5.1)遍历选择一个子高斯分布作为当前的子高斯分布i；

5.2)对子高斯分布i，根据式(14)将初始状态相对标称值的初始状态偏差δx0变换成子高斯分布i的均值μi和相应偏差δxi；

δx0＝μi+δxi(14)

式(14)中，δx0为初始状态偏差，μi为子高斯分布i的均值，δxi为子高斯分布i的偏差；

5.3)将变换成子高斯分布i的均值μi和相应偏差δxi的初始状态偏差δx0代入终端状态关于初始状态偏差δx0的高阶泰勒展开多项式中，得到式(15)所示终端状态；

式(15)中，xf为终端状态，t为高阶泰勒展开多项式，t的上标n为偏差演化阶数、t的下标xf为终端状态，δxi为子高斯分布i的偏差；

5.4)根据式(16)计算终端时刻子高斯分布i的均值μi_f作为式(15)所示终端状态的期望值；

式(16)中，μi_f为终端时刻子高斯分布i的均值，xf为终端状态，t为高阶泰勒展开多项式，t的上标n为偏差演化阶数、t的下标xf为终端状态，δxi为子高斯分布i的偏差；

5.5)根据式(17)将终端时刻子高斯分布i的协方差矩阵pi_f表示成关于初始子高斯分布偏差δxi的2n阶泰勒展开式并求出子高斯分布i对应的期望值；

式(17)中，pi_f表示终端时刻子高斯分布i的协方差矩阵，xf为终端状态，μi_f为终端时刻子高斯分布i的均值，t为高阶泰勒展开多项式，t的上标中的n为偏差演化阶数，δxi为子高斯分布i的偏差；

5.6)判断所有的子高斯分布是否已经遍历完毕，如果尚未遍历完毕，则跳转执行步骤5.2)；否则，跳转执行下一步；

5.7)在得到所有子高斯分布在终端时刻的均值μi_f和协方差矩阵pi_f的基础上，根据式(18)求取终端时刻系统状态分布的精确概率密度函数并输出；

式(18)中，p(t,xf)表示终端时刻系统状态分布的精确概率密度函数，t为终端时刻，xf为终端状态，wi为第i个子高斯分布对应的权重，m为打靶次数，pg(xf；μi_f,pi_f)为第i个子高斯分布的概率密度函数，μi_f为终端时刻子高斯分布i的均值，pi_f表示终端时刻子高斯分布i的协方差矩阵。

本发明基于微分代数与高斯和的非线性系统状态偏差演化方法具有下述优点：

1、本发明基于微分代数方法，能够以任意指定阶精度预报非线性系统的状态，而且采用了一种递推算法计算任意阶的高斯分布中心矩，因此基于微分代数与高斯和的非线性系统状态偏差演化方法可以自动拓展至任意指定阶偏差演化精度，无需手动推导动力学方程的高阶偏导数，具有使用方便、计算精度高的优点，适用于非线性强和长时间预报的非线性系统偏差演化分析问题。

2、本发明采用高斯和模型表征系统状态偏差分布的概率密度函数，具有对非线性系统偏差非高斯性的精确描述能力，而且各子高斯分布预报均基于同一终端状态关于初始状态偏差的高阶泰勒展开多项式，不需为每个子高斯分布分别预报终端状态，因此基于微分代数与高斯和的非线性系统状态偏差演化方法在精确描述偏差分布及其非高斯性的同时，仍能保持相对montecarlo仿真方法的显著效率优势。

3、在航天器轨道动力学领域，本发明可应用于航天器绝对或相对轨迹预报及其偏差演化分析、轨道确定误差分析等，获得的偏差演化信息可进一步用于空间碰撞预警、航天器轨迹安全性分析与设计等，具有使用方便、计算精度和计算效率高等优点，适用于强非线性系统的长时间偏差演化分析问题。

附图说明

图1为本发明实施例一方法的基本流程示意图。

图2为本发明实施例一中求子高斯分布协方差矩阵pi的流程示意图。

图3为本发明实施例一的终端状态偏差分布对比图。

图4为本发明实施例一的终端概率密度函数对比图。

图5为本发明实施例二的终端状态偏差分布对比图。

图6为本发明实施例二的终端概率密度函数对比图。

具体实施方式

实施例一：

如图1所示，本实施例基于微分代数与高斯和的非线性系统状态偏差演化方法步骤包括：

1)输入非线性系统的动力学方程、偏差演化阶数n、初始状态相对标称值和初始偏差协方差矩阵p0；

2)采用微分代数方法，针对非线性系统的动力学方程预报非线性系统的终端状态，并表示为关于初始状态偏差δx0的高阶泰勒展开多项式；

3)针对初始偏差协方差矩阵p0确定子高斯分布协方差矩阵pi，针对子高斯分布协方差矩阵pi中的多个子高斯分布进行打靶拟合初始偏差分布生成高斯和模型；

4)针对生成高斯和模型计算给定阶数的子高斯分布高阶中心矩；

5)基于终端状态关于初始状态偏差的高阶多项式和子高斯分布的各阶中心矩，确定各子高斯分布在终端时刻的均值和协方差矩阵，进而给出高斯和形式的终端状态偏差分布概率密度函数。

本实施例基于微分代数与高斯和的非线性系统状态偏差演化方法采用微分代数方法，预报非线性系统终端状态，并将其表示为关于初始状态偏差的高阶多项式；采用高斯和模型，以多个子高斯分布加权和的形式拟合初始偏差分布；利用子高斯分布的协方差矩阵，计算给定阶数的子高斯分布高阶中心矩；基于终端状态关于初始状态偏差的高阶多项式和子高斯分布的各阶矩，确定各子高斯分布在终端时刻的均值和协方差矩阵。本发明面向一般的非线性动力学系统，动力学模型可以是状态转移方程形式，或常微分方程形式，可自动进行任意高阶的状态预报及其偏差演化分析。在航天器轨道动力学领域，本发明可具体应用于航天器绝对或相对轨迹预报及其偏差演化分析、轨道确定误差分析等，获得的偏差演化信息可进一步用于空间碰撞预警、航天器轨迹安全性分析与设计等，具有使用方便、计算精度和计算效率高等优点，适用于强非线性系统的长时间偏差演化分析问题。

本实施例基于微分代数与高斯和的非线性系统状态偏差演化方法面向一般的非线性动力学系统，动力学模型可以是状态转移方程形式，或常微分方程形式，可自动进行任意高阶的状态预报及其偏差演化分析。本实施例中，步骤1)中非线性系统的动力学方程为非线性映射形式xf＝f(x0)或常微分方程形式且所述非线性映射形式的函数f(·)或者常微分方程形式的函数g(·)均为任意从n维实数域rⁿ到rⁿ的非线性映射。

本实施例中非线性系统为近地航天器开普勒轨道，系统维数n＝6，状态变量x为地球惯性坐标系下航天器的位置和速度矢量，即x＝[r^tv^t]^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可以计算出初始时刻航天器标称状态给定轨道预报时间tf为一个轨道周期，即tf＝5553.63s。

此外，本实施例给定偏差演化阶数n＝4，给定的初始偏差协方差矩阵p0为：

本实施例中，步骤2)的详细步骤包括：

2.1)将初始状态x0初始化为式(1)所示包含初始状态相对标称值的偏差形式；

式(1)中，x0为初始状态，为初始状态相对标称值，δx0为初始偏差；

2.2)将初始状态x0代入非线性系统的动力学方程，应用微分代数方法，求而出关于初始偏差δx0的高阶泰勒展开多项式描述的终端状态如式(2)所示；

式(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)所示函数表达式计算矩阵的奇异值分解，其中r0为初始偏差协方差矩阵，rmin为子高斯分布的平方根信息矩阵的下界；

式(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；

式(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为子高斯分布的平方根信息矩阵的下界；

式(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)所示；

式(8)中，δyi为协方差减小量的平方根，r0为初始偏差协方差矩阵，δri为信息变化矩阵，矩阵rdi通过如式(9)所示qr分解得到；

式(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)所示；

式(10)中，p(t,x0)为高斯和形式的初始偏差概率密度函数，wi为第i个子高斯分布对应的权重，m为打靶次数，pg(x0；μi,pi)为第i个子高斯分布的概率密度函数，且第i个子高斯分布的概率密度函数pg(x0；μi,pi)的函数表达式如式(11)所示；

式(11)中，n为系统维数，pi为第i个子高斯分布对应的协方差矩阵，x0为初始状态，μi为第i个子高斯分布的均值。

本实施例中，步骤4)的详细步骤包括：

4.1)根据高斯分布的特性，确定子高斯分布的任意奇数阶中心矩均为0；

例如，子高斯分布的一阶中心矩为：

e{δx^a}＝0

其中，e{·}表示求括号内物理量的期望值，e{δx^a}为子高斯分布偏差的一阶中心矩，x^a是状态变量x的第a个元素，δx^a为状态变量x与状态均值μi第a个元素的偏差。

例如，子高斯分布的三阶中心矩为：

e{δx^aδx^bδx^c}＝0,a,b,c∈[1,2,...,n]

其中，e{·}表示求括号内物理量的期望值，e{δx^aδx^bδx^c}为子高斯分布偏差的一阶中心矩，x^a是状态变量x的第a个元素，δx^a为状态变量x与状态均值μi第a个元素的偏差，x^b为状态变量x的第b个元素，δx^b为状态变量x与状态均值μi第b个元素的偏差，x^c为状态变量x的第c个元素，δx^c为状态变量x与状态均值μi第c个元素的偏差。更高奇数阶中心矩类似，均为0。

4.2)对第i个子高斯分布而言，令其中a＝1,2,...,n，x^a是状态变量x的第a个元素，是状态均值μi的第a个元素；根据式(12)计算子高斯分布偏差的二阶中心矩；

e{δx^aδx^b}＝pi^ab,a,b∈[1,2,...,n](12)

式(12)中，e{·}表示求括号内物理量的期望值，e{δx^aδx^b}为子高斯分布偏差的二阶中心矩，δx^a为状态变量x与状态均值μi第a个元素的偏差，δx^b为状态变量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为偏差演化阶数；

式(13)中，e{·}表示求括号内物理量的期望值，e{δx^aδx^bδx^c...δx^d}为子高斯分布偏差的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{δx^c...δx^d}为子高斯分布偏差的m阶中心矩的第(c,…,d)个元素；pi^ac表示第i个子高斯分布的协方差矩阵pi的第a行第c列元素，e{δx^b...δx^d}为子高斯分布偏差的m阶中心矩的第(b,…,d)个元素；以此类推，直至pi^ad表示第i个子高斯分布的协方差矩阵pi的第a行第d列元素，e{δx^bδx^c...}为子高斯分布偏差的m阶中心矩的第(b,c,…)个元素。

本实施例中，步骤5)的详细步骤包括：

5.1)遍历选择一个子高斯分布作为当前的子高斯分布i；

5.2)对子高斯分布i，根据式(14)将初始状态相对标称值的初始状态偏差δx0变换成子高斯分布i的均值μi和相应偏差δxi；

δx0＝μi+δxi(14)

式(14)中，δx0为初始状态偏差，μi为子高斯分布i的均值，δxi为子高斯分布i的偏差；

5.3)将变换成子高斯分布i的均值μi和相应偏差δxi的初始状态偏差δx0代入终端状态关于初始状态偏差δx0的高阶泰勒展开多项式中，得到式(15)所示终端状态；

式(15)中，xf为终端状态，t为高阶泰勒展开多项式，t的上标n为偏差演化阶数、t的下标xf为终端状态，δxi为子高斯分布i的偏差；

5.4)根据式(16)计算终端时刻子高斯分布i的均值μi_f作为式(15)所示终端状态的期望值；

式(16)中，μi_f为终端时刻子高斯分布i的均值，xf为终端状态，t为高阶泰勒展开多项式，t的上标n为偏差演化阶数、t的下标xf为终端状态，δxi为子高斯分布i的偏差；

5.5)根据式(17)将终端时刻子高斯分布i的协方差矩阵pi_f表示成关于初始子高斯分布偏差δxi的2n阶泰勒展开式并求出子高斯分布i对应的期望值；

式(17)中，pi_f表示终端时刻子高斯分布i的协方差矩阵，xf为终端状态，μi_f为终端时刻子高斯分布i的均值，t为高阶泰勒展开多项式，t的上标中的n为偏差演化阶数，δxi为子高斯分布i的偏差；

5.6)判断所有的子高斯分布是否已经遍历完毕，如果尚未遍历完毕，则跳转执行步骤5.2)；否则，跳转执行下一步；

5.7)在得到所有子高斯分布在终端时刻的均值μi_f和协方差矩阵pi_f的基础上，根据式(18)求取终端时刻系统状态分布的精确概率密度函数并输出；

式(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＝[r^tv^t]^t，其中r＝[x,y,z]^t是相对位置矢量，v＝[vx,vy,vz]^t是相对速度矢量。

本实施例中，非线性系统的动力学方程为常微分方程形式，终端状态通过变步长龙格库塔积分方法得到。

本实施例中，非线性系统的动力学方程的具体表达式为：

其中，a为主航天器轨道半长轴，本实施例中即为地球静止轨道半径a＝42164200m，静止轨道航天器的轨道角速度为n＝7.2921×10^-5rad/s，地球引力常数为μ＝3.986×10¹⁴m³/s²，太阳光压摄动加速度为γ＝[γx,γy,γz]^t，其具体表达形式可参考文献：fehse,w.,“disturbancebysolarpressureofrendezvoustrajectoriesingeo,”5thicatt,noordwijk,netherlands,may2012。

给定初始时刻航天器经典轨道根数标称值为：

x0＝[0,-70000,0,0,-1,0](m,m/s)

给定轨道预报时间tf为一个轨道周期，本实施例中即为tf＝86164s。

步骤2.1)中，给定初始状态协方差矩阵p0为：

相比实施例一而言，本实施例非线性系统变为地球静止轨道考虑太阳光压的航天器非线性相对运动轨迹，动力学方程变为常微分方程形式，终端状态通过变步长龙格库塔积分方法得到。

参见图5所示终端状态偏差分布对比图中，子高斯分布表示本发明预报得到的终端时刻子高斯偏差分布情况，montecarlo打靶点表示传统montecarlo仿真方法得到的终端时刻打靶点分布情况，montecarlo协方差表示根据montecarlo打靶点统计出的协方差3-sigma椭球，montecarlo均值表示根据montecarlo打靶点统计出的分布均值，线性化协方差表示传统线性化方法得到的终端时刻协方差3-sigma椭球，线性化均值表示传统线性化方法得到的终端时刻分布均值。参见图6所示终端概率密度函数对比图中，高斯和表示本发明得到的终端高斯和模型概率密度函数，montecarlo表示根据montecarlo打靶点统计出的概率密度函数，线性化方法表示传统线性化方法得到的终端时刻概率密度函数。参见图5所示终端状态偏差分布对比和图6所示终端概率密度函数对比，本实施例基于微分代数与高斯和的非线性系统状态偏差演化方法仍然能很好地拟合了montecarlo仿真结果，有效捕获航天器相对轨迹偏差分布的非高斯性，并准确预报终端时刻偏差分布的概率密度函数。本实施例验证了本实施例基于微分代数与高斯和的非线性系统状态偏差演化方法的普适性，对一般的非线性系统，动力学方程可以是非线性映射形式或常微分方程形式，均能精确预报其状态偏差的演化情况。

以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。