CN106970643B - 一种解析的卫星非线性相对运动偏差传播分析方法 - Google Patents

一种解析的卫星非线性相对运动偏差传播分析方法 Download PDF

Info

Publication number
CN106970643B
CN106970643B CN201710287256.9A CN201710287256A CN106970643B CN 106970643 B CN106970643 B CN 106970643B CN 201710287256 A CN201710287256 A CN 201710287256A CN 106970643 B CN106970643 B CN 106970643B
Authority
CN
China
Prior art keywords
satellite
relative motion
deviation
state
time
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
CN201710287256.9A
Other languages
English (en)
Other versions
CN106970643A (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 CN201710287256.9A priority Critical patent/CN106970643B/zh
Publication of CN106970643A publication Critical patent/CN106970643A/zh
Application granted granted Critical
Publication of CN106970643B publication Critical patent/CN106970643B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05D—SYSTEMS FOR CONTROLLING OR REGULATING NON-ELECTRIC VARIABLES
    • G05D1/00—Control of position, course, altitude or attitude of land, water, air or space vehicles, e.g. using automatic pilots
    • G05D1/10—Simultaneous control of position or course in three dimensions
    • G05D1/101—Simultaneous control of position or course in three dimensions specially adapted for aircraft
    • G05D1/104—Simultaneous control of position or course in three dimensions specially adapted for aircraft involving a plurality of aircrafts, e.g. formation flying

Landscapes

  • Engineering & Computer Science (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明公开了一种解析的卫星非线性相对运动偏差传播分析方法，步骤包括指定主卫星与从卫星，并输入初始时刻的参考卫星绝对轨道状态、两卫星标称相对运动状态、初始相对运动状态偏差的概率密度函数，根据考虑J2摄动的非线性相对运动方程计算用于解析预报卫星相对运动状态及偏差的一阶及二阶状态转移张量，如为高斯分布则采用协方差分析方法计算并输出解析结果，否则采用高斯和模型计算并输出解析结果。本发明考虑了J2摄动项和二阶非线性项，能够用于两相距较远卫星的相对运动状态偏差的长时间、高精度解析预报，获得的偏差信息可用于编队卫星碰撞概率计算及碰撞预警，具有设计方法正确合理、方法解析计算量小、对实际工程任务的适用性好的优点。

Description

一种解析的卫星非线性相对运动偏差传播分析方法
技术领域
本发明涉及航天器相对运动状态偏差传播分析技术，具体涉及一种解析的卫星非线性相对运动偏差传播分析方法。
背景技术
航天器轨道偏差传播分析在空间态势感知各功能任务(如目标追踪与数据关联、碰撞预警、传感器导引、机动检测等)中具有重要应用。随着空间目标不断增加，而观测设备有限，大多数空间目标不能连续跟踪，因此需要对其轨道及偏差进行长时间预报，以引导下一次跟踪或对可能的碰撞进行预警。
近年来，卫星编队飞行在理论研究和工程实践中都受到广泛关注，因为通过多个航天器编队飞行，可以低成本、高效率的完成很多航天任务，如构成合成孔径雷达、太空高精度干涉测量、对地观测等。编队卫星的构型设计及保持需要知道卫星的绝对运动状态及卫星间的相对运动状态，由于导航设备误差、星间通信链路时延、及观测盲区等因素影响，使获得的编队卫星的相对运动状态存在不确定性。由于编队卫星间距离相对较近，如果不对存在于相对运动状态中的偏差进行精确监控，极有可能引起编队卫星相互碰撞。因此，偏差传播分析结果作为碰撞概率计算的重要依据，对编队、集群飞行卫星的碰撞预警及规避任务尤为重要。
轨道动力学中经典的偏差演化分析方法有线性协方差分析方法和Monte Carlo仿真方法。线性协方差分析是通过将动力系统线性化的方式来进行偏差传播，该方法使用简单、计算量小，但是对强非线性系统或长时间偏差传播问题计算误差较大。Monte Carlo仿真虽然能通过统计的方法得到轨道偏差的高精度分布属性，但需要进行大量抽样仿真，计算量大。近年来，不少学者提出了一些高精度半解析或数值的偏差演化分析方法，有效兼顾了偏差演化分析中的计算精度与计算效率。然而，这些方法多是针对卫星绝对状态偏差的传播，且需要进行一定量的轨道积分、依然有较大的计算量，对卫星编队飞行任务不太适用。
特别地，对编队飞行卫星的偏差分析与碰撞预警多是在轨进行的。由于星载计算机的计算能力相对较弱，因此需要偏差传播分析方法能解析计算，以便星载计算机能快速、准确地进行碰撞预警。因此，如何基于解析的卫星非线性相对运动偏差传播分析，以高精度地对卫星相对运动偏差进行长时间预报，已经成为一项亟待解决的关键技术问题。
发明内容
本发明要解决的技术问题：针对现有技术的上述问题，提供一种基于考虑J2摄动的非线性相对运动方程，能够用于两相距较远卫星的相对运动状态偏差的长时间、高精度解析预报，获得的偏差信息可用于编队卫星碰撞概率计算及碰撞预警，设计方法正确合理、对实际工程任务的适用性好的解析的卫星非线性相对运动偏差传播分析方法。
为了解决上述技术问题，本发明采用的技术方案为：
一种解析的卫星非线性相对运动偏差传播分析方法，实施步骤包括：
1)为编队飞行的两卫星指定主卫星与从卫星，并输入初始时刻的参考卫星绝对轨道状态、两卫星标称相对运动状态、输入初始相对运动状态偏差的概率密度函数；
2)基于初始时刻的参考卫星绝对轨道状态、两卫星标称相对运动状态，根据考虑J2摄动的非线性相对运动方程计算用于解析预报卫星相对运动状态的一阶及二阶状态转移张量；
3)基于卫星相对运动状态的一阶及二阶状态转移张量，计算用于解析预报卫星相对运动状态偏差的一阶及二阶状态转移张量；
4)判断初始及终端时刻的卫星偏差是否均为高斯分布，如果均为高斯分布则跳转执行步骤5)，否则，跳转执行步骤6)；
5)采用协方差分析方法，获得卫星相对运动状态偏差的均值及协方差矩阵的解析非线性传播方程，跳转执行步骤7)；
6)采用高斯和模型，按步骤6.1)～6.3)获得卫星非线性相对运动状态偏差的概率密度函数的解析传播方程，然后跳转执行步骤7)；
6.1)用N个子高斯分布概率密度函数去逼近初始相对状态偏差的概率密度函数p(t0,dx0)获得每一个子高斯分布的权重ωi、均值mi(t0)和协方差矩阵Ci(t0)；
6.2)采用卫星相对运动状态偏差的均值及协方差矩阵的解析非线性传播方程，将每个子高斯分布的均值mi(t0)及协方差矩阵Ci(t0)预报到终端时刻tf，预报过程中子高斯分布的权重ωi保持不变，获得卫星终端时刻每个子高斯的均值mi(tf)及协方差矩阵Ci(tf)，其中i＝1,2,…,N，N为子高斯分布概率密度函数的数量；
6.3)针对第i个子高斯分布，其中i＝1,2,…,N，N为子高斯分布概率密度函数的数量，根据其权重系数ωi、均值mi(tf)及协方差矩阵Ci(tf)，获得卫星非线性相对运动状态偏差的解析传播均值为
Figure GDA0002222800330000021
协方差矩阵为
Figure GDA0002222800330000022
概率密度函数为
Figure GDA0002222800330000023
其中
Figure GDA0002222800330000024
dxf为卫星终端相对运动状态偏差，pg(dxf；mi(tf),Ci(tf))表示tf时刻第i个子高斯分布的概率密度函数；
7)输出终端时刻两卫星相对运动状态偏差的均值、协方差矩阵及概率密度函数。
优选地，步骤1)中输入初始时刻的参考卫星绝对轨道状态具体是指输入主卫星初始时刻t0的绝对轨道状态E0＝[a0,e0,i0,Ω0,ω0,f0]，其中，E0为主卫星初始时刻t0的绝对轨道状态，a0为主卫星初始时刻t0的半长轴，e0为主卫星初始时刻t0的偏心率，i0为主卫星初始时刻t0的轨道倾角，Ω0为主卫星初始时刻t0的升交点赤经，ω0为主卫星初始时刻t0的近拱点角距，f0为主卫星初始时刻t0的真近点角；且输入初始时刻t0从卫星相对主卫星的标称相对运动状态为δx(t0)，δx(t0)在主卫星当地轨道坐标系中表示，该坐标系原点为主卫星质心，x轴沿主卫星地心矢径方向，z轴沿轨道面法向，y轴与x、z轴构成右手坐标系；输入初始相对运动状态偏差为dx(t0)及其概率密度函数为p(t0,dx0)，若初始偏差为高斯分布，则其概率密度函数为pg(dx0；m0,C0)，m0为初始相对状态偏差的均值矩阵，C0为初始相对状态偏差的协方差矩阵；
优选地，步骤2)中考虑J2摄动的非线性相对运动方程的函数表达式如式(1)所示；
Figure GDA0002222800330000031
式(1)中，δx(tf)表示时刻tf从卫星相对主卫星的标称相对运动状态，Φ(tf,t0)是从t0时刻到tf时刻的一阶状态转移张量，δx(t0)表示初始时刻t0从卫星相对主卫星的标称相对运动状态，Ψ(tf,t0)为从t0时刻到tf时刻的二阶状态转移张量，T(tf)为以纬度幅角θf为自变量的无量纲化坐标到以时间tf为自变量的量纲化坐标的转换矩阵，
Figure GDA0002222800330000032
为从θ0纬度幅角到θf纬度幅角的状态转移矩阵，T(t0)为以纬度幅角θ0为自变量的无量纲化坐标到以时间t0为自变量的量纲化坐标的转换矩阵，
Figure GDA0002222800330000033
为从θ0纬度幅角到θf纬度幅角的二阶状态转移张量，T-1(t0)为转换矩阵T(t0)的逆矩阵，
Figure GDA0002222800330000034
为Kronecker张量积；
优选地，步骤2)中计算用于解析预报卫星相对运动状态的一阶状态转移张量如式(2)所示，计算用于解析预报卫星相对运动状态的二阶状态转移张量如式(3)所示；
Figure GDA0002222800330000035
Figure GDA0002222800330000036
式(2)和式(3)中，
Figure GDA0002222800330000043
为用于解析预报卫星相对运动状态的从θ0纬度幅角到θf纬度幅角的一阶状态转移张量，
Figure GDA0002222800330000044
表示θf＝θ0时的中间变量
Figure GDA0002222800330000045
Figure GDA0002222800330000046
表示
Figure GDA0002222800330000047
的逆矩阵，
Figure GDA0002222800330000048
表示θf＝θ0时的中间变量
Figure GDA0002222800330000049
Figure GDA00022228003300000410
为用于解析预报卫星相对运动状态的从θ0纬度幅角到θf纬度幅角的二阶状态转移张量，
Figure GDA00022228003300000411
为Kronecker张量积，
Figure GDA00022228003300000412
和
Figure GDA00022228003300000413
均为公式推导使用的中间变量，且中间变量
Figure GDA00022228003300000414
和
Figure GDA00022228003300000415
的计算函数式如式(4)所示；
Figure GDA0002222800330000041
式(4)中，A、B为公式推导使用的中间变量，
Figure GDA00022228003300000416
表示t0时刻平均轨道要素偏差到tf时刻平均轨道要素偏差的一阶传递矩阵，
Figure GDA00022228003300000417
表示t0时刻轨道要素偏差的无量纲化矩阵，
Figure GDA00022228003300000418
表示tf时刻轨道要素偏差的无量纲化矩阵，
Figure GDA00022228003300000419
表示t0时刻平均轨道要素偏差到tf时刻平均轨道要素偏差的二阶传递张量，Q(θf)表示tf时刻无量纲轨道要素偏差到相对状态的二阶转换张量，T(θf)表示tf时刻以纬度幅角为自变量的无量纲化坐标到以时间为自变量的量纲化坐标的转换矩阵，T-1(θf)表示T(θf)的逆矩阵，Π表示相对状态分量次序变换矩阵，Σ(θf)表示tf时刻轨道要素偏差到相对状态的转换矩阵，
Figure GDA00022228003300000420
表示平均轨道要素偏差到吻切轨道要素偏差的转换矩阵，
Figure GDA00022228003300000421
为Kronecker张量积。
优选地，步骤3)计算用于解析预报卫星相对运动状态偏差的一阶及二阶状态转移张量的函数表达式如式(5)所示；
Figure GDA0002222800330000042
式(5)中，dxf为终端时刻的相对运动状态偏差，
Figure GDA00022228003300000422
为用于解析预报卫星相对运动状态偏差的一阶状态转移张量，dx0为初始的相对运动状态偏差，
Figure GDA00022228003300000423
为用于解析预报卫星相对运动状态偏差的二阶状态转移张量，
Figure GDA00022228003300000424
为
Figure GDA00022228003300000425
的第i行第j列元素，
Figure GDA00022228003300000426
为
Figure GDA00022228003300000427
的第i行第j列元素，
Figure GDA00022228003300000428
为
Figure GDA00022228003300000429
的第i行k列j维元素，
Figure GDA00022228003300000430
是从t0时刻到tf时刻的状态转移矩阵，
Figure GDA0002222800330000053
为从t0时刻到tf时刻的二阶状态转移张量，
Figure GDA0002222800330000054
为初始时刻t0从卫星相对主卫星的标称相对运动状态δx0的第k个分量，
Figure GDA0002222800330000055
为Kronecker张量积。
优选地，步骤5)中获得的卫星相对运动状态偏差的均值及协方差矩阵的解析非线性传播方程的函数表达式如式(6)所示；
Figure GDA0002222800330000051
式(6)中，mi(tf)为卫星相对状态偏差的均值矩阵m(tf)的第i行元素，mj(tf)为卫星相对状态偏差的均值矩阵m(tf)的第j行元素，Cij(tf)为卫星相对状态偏差的矩阵，
Figure GDA0002222800330000056
为用于解析预报卫星相对运动状态偏差的一阶状态转移张量
Figure GDA0002222800330000057
的第i行第a列元素，
Figure GDA0002222800330000058
为用于解析预报卫星相对运动状态偏差的一阶状态转移张量
Figure GDA0002222800330000059
的第j行第b列元素，
Figure GDA00022228003300000510
为从t0时刻到tf时刻的二阶状态转移张量
Figure GDA00022228003300000511
的第i行a列b维元素，
Figure GDA00022228003300000512
为从t0时刻到tf时刻的二阶状态转移张量
Figure GDA00022228003300000513
的第j行b列c维元素，
Figure GDA00022228003300000514
为从t0时刻到tf时刻的二阶状态转移张量
Figure GDA00022228003300000515
的第i行b列c维元素，
Figure GDA00022228003300000516
为从t0时刻到tf时刻的二阶状态转移张量
Figure GDA00022228003300000517
的第j行c列d维元素，
Figure GDA00022228003300000518
分别为初始相对状态偏差dx0的前四阶矩。
优选地，初始相对状态偏差dx0的前四阶矩的计算函数表达式如式(7)所示；
Figure GDA0002222800330000052
式(7)中，
Figure GDA00022228003300000519
分别为初始相对状态偏差dx0的前四阶矩，
Figure GDA00022228003300000520
分别为初始相对状态偏差的均值矩阵m0的a,b,c,d分量元素，
Figure GDA00022228003300000521
分别为初始相对状态偏差的协方差矩阵C0的ab分量元素、bc分量元素、cd分量元素、ad分量元素。
优选地，步骤6)的详细步骤包括：
6.1)用N个子高斯分布概率密度函数去逼近初始相对状态偏差的概率密度函数p(t0,dx0)获得每一个子高斯分布的权重、均值和协方差矩阵；
6.2)采用卫星相对运动状态偏差的均值及协方差矩阵的解析非线性传播方程，将每个子高斯分布的均值mi(t0)及协方差矩阵Ci(t0)预报到终端时刻tf，预报过程中子高斯分布的权重保持不变，获得终端时刻每个子高斯的均值mi(tf)及协方差矩阵Ci(tf)，其中i＝1,2,…,N，N为高斯分布概率密度函数的数量；
6.3)针对每一个第i个子高斯分布，其中i＝1,2,…,N，N为子高斯分布概率密度函数的数量，根据其权重系数ωi、均值mi(tf)及协方差矩阵Ci(tf)，根据式(8)和式(9)所示获得卫星非线性相对运动状态偏差的概率密度函数的解析传播方程计算卫星相对状态偏差的均值及协方差矩阵以及卫星偏差的概率密度函数p(tf,dxf)；
Figure GDA0002222800330000061
Figure GDA0002222800330000062
式(8)和(9)中，m(tf),C(tf)为卫星相对状态偏差dxf的均值及协方差矩阵，ωi表示第i个子高斯分布的权重系数，N为高斯分布概率密度函数的数量，ω为所有高斯分布的权重系数之和，mi(tf)为tf时刻第i个子高斯分布的均值，Ci(tf)为tf时刻第i个子高斯分布的协方差矩阵，
Figure GDA0002222800330000064
为mi(tf)的转置矩阵，mT(tf)为m(tf)的转置矩阵，p(tf,dxf)为卫星偏差的概率密度函数，pg(dxf；mi(tf),Ci(tf))表示tf时刻第i个子高斯分布的概率密度函数。
优选地，步骤6.1)中逼近初始相对状态偏差的概率密度函数p(t0,dx0)采用的逼近公式如式(10)所示；
Figure GDA0002222800330000063
式(10)中，p(t0,dx0)为初始相对状态偏差的概率密度函数，dx0为初始相对状态偏差，pg(dx0；mi,Ci)表示第i个子高斯分布的概率密度函数，mi和Ci分别为相应的均值和协方差矩阵，ωi表示第i个子高斯分布的权重系数。
优选地，第i个子高斯分布的权重系数ωi的计算函数表达式如式(11)所示；
Figure GDA0002222800330000071
式(11)中，ωi表示第i个子高斯分布的权重系数，m0、C0为初始相对状态偏差的均值及协方差矩阵，
Figure GDA0002222800330000072
及
Figure GDA0002222800330000073
为将一维标准正态分布pg(x；0,1)等方差地分割为N个子正态分布的权重、均值及标准差，diag{a}表示以a向量为对角线元素的矩阵，V和Λ分别为对C0进行特征值分解所得的特征向量矩阵与特征值矩阵，λ1,λ2,...,λn为C0的特征值，λj与υj分别为沿选定的j方向的特征值与特征向量，即Λ的第j行j列元素和V的第j列向量，j选择为初始偏差协方差矩阵中径向位置或横向速度的对应值。
本发明解析的卫星非线性相对运动偏差传播分析方法具有下述优点：1、本发明通过计算相对状态偏差的状态转移矩阵及张量，计算相对运动状态偏差均值、协方差矩阵、及概率密度函数的非线性传播结果，考虑了航天器实际飞行环境中的主要摄动因素J2项及二阶非线性项，可用于两相距较远航天器的长时间、高精度相对运动状态偏差传播分析。2、本发明利用编队卫星的初始相对运动状态均值及协方差矩阵，计算相对状态及相对状态偏差的状态转移矩阵及张量，预报相对运动状态偏差的均值、协方差矩阵、及概率密度函数，方法解析，计算速度快，可用于卫星在轨碰撞预警及碰撞规避。
附图说明
图1为本发明实施例一方法的基本流程示意图。
图2为本发明实施例一方法和和精确的蒙特卡洛仿真的终端误差分布对比图。
图3为本发明实施例一方法用的高斯和模型对均值和协方差的预报与蒙特卡洛仿真的样本点分布对比图。
图4为本发明实施例二中的终端误差均值及协方差矩阵的预报结果对比图。
具体实施方式
实施例一：
如图1所示，本实施例解析的卫星非线性相对运动偏差传播分析方法的实施步骤包括：
1)为编队飞行的两卫星指定主卫星与从卫星，并输入初始时刻的参考卫星绝对轨道状态、两卫星标称相对运动状态、输入初始相对运动状态偏差的概率密度函数；
2)基于初始时刻的参考卫星绝对轨道状态、两卫星标称相对运动状态，根据考虑J2摄动的非线性相对运动方程计算用于解析预报卫星相对运动状态的一阶及二阶状态转移张量；
3)基于卫星相对运动状态的一阶及二阶状态转移张量，计算用于解析预报卫星相对运动状态偏差的一阶及二阶状态转移张量；
4)判断初始及终端时刻的卫星偏差是否均为高斯分布，如果均为高斯分布则跳转执行步骤5)，否则，跳转执行步骤6)；
5)采用协方差分析方法，获得卫星相对运动状态偏差的均值及协方差矩阵的解析非线性传播方程，跳转执行步骤7)；
6)采用高斯和模型，获得卫星非线性相对运动状态偏差的概率密度函数的解析传播方程，跳转执行步骤7)；
7)输出终端时刻两卫星相对运动状态偏差的均值、协方差矩阵及概率密度函数。
本实施例解析的卫星非线性相对运动偏差传播分析方法通过基于考虑J2摄动的非线性相对运动方程，能够用于两相距较远卫星的相对运动状态偏差的长时间、高精度解析预报，获得的偏差信息可用于编队卫星碰撞概率计算及碰撞预警，具有设计方法正确合理、方法解析计算量小、对实际工程任务的适用性好的优点。
本实施例中，步骤1)中输入初始时刻的参考卫星绝对轨道状态具体是指输入主卫星初始时刻t0的绝对轨道状态E0＝[a0,e0,i0,Ω0,ω0,f0]，其中，E0为主卫星初始时刻t0的绝对轨道状态，a0为主卫星初始时刻t0的半长轴，e0为主卫星初始时刻t0的偏心率，i0为主卫星初始时刻t0的轨道倾角，Ω0为主卫星初始时刻t0的升交点赤经，ω0为主卫星初始时刻t0的近拱点角距，f0为主卫星初始时刻t0的真近点角；且输入初始时刻t0从卫星相对主卫星的标称相对运动状态为δx(t0)，δx(t0)在主卫星当地轨道坐标系中表示，该坐标系原点为主卫星质心，x轴沿主卫星地心矢径方向，z轴沿轨道面法向，y轴与x、z轴构成右手坐标系；输入初始相对运动状态偏差为dx(t0)及其概率密度函数为p(t0,dx0)，若初始偏差为高斯分布，则其概率密度函数为pg(dx0；m0,C0)，m0为初始相对状态偏差的均值矩阵，C0为初始相对状态偏差的协方差矩阵；
本实施例中，步骤2)中考虑J2摄动的非线性相对运动方程的函数表达式如式(1)所示；
Figure GDA0002222800330000081
式(1)中，δx(tf)表示时刻tf从卫星相对主卫星的标称相对运动状态，Φ(tf,t0)是从t0时刻到tf时刻的状态转移矩阵，δx(t0)表示初始时刻t0从卫星相对主卫星的标称相对运动状态，Ψ(tf,t0)为从t0时刻到tf时刻的二阶状态转移张量，T(tf)为以纬度幅角θf为自变量的无量纲化坐标到以时间tf为自变量的量纲化坐标的转换矩阵，
Figure GDA0002222800330000094
为从θ0纬度幅角到θf纬度幅角的一阶状态转移张量，T(t0)为以纬度幅角θ0为自变量的无量纲化坐标到以时间t0为自变量的量纲化坐标的转换矩阵，
Figure GDA0002222800330000095
为从θ0纬度幅角到θf纬度幅角的二阶状态转移张量，T-1(t0)为转换矩阵T(t0)的逆矩阵，
Figure GDA0002222800330000096
为Kronecker张量积；
本实施例中，步骤2)中计算用于解析预报卫星相对运动状态的一阶状态转移张量如式(2)所示，计算用于解析预报卫星相对运动状态的二阶状态转移张量如式(3)所示；
Figure GDA0002222800330000091
Figure GDA0002222800330000092
式(2)和式(3)中，
Figure GDA0002222800330000097
为用于解析预报卫星相对运动状态的从θ0纬度幅角到θf纬度幅角的一阶状态转移张量，
Figure GDA0002222800330000098
表示θf＝θ0时的中间变量
Figure GDA0002222800330000099
Figure GDA00022228003300000910
表示
Figure GDA00022228003300000911
的逆矩阵，
Figure GDA00022228003300000917
表示θf＝θ0时的中间变量
Figure GDA00022228003300000912
Figure GDA00022228003300000913
为用于解析预报卫星相对运动状态的从θ0纬度幅角到θf纬度幅角的二阶状态转移张量，
Figure GDA00022228003300000914
为Kronecker张量积，
Figure GDA00022228003300000915
和
Figure GDA00022228003300000916
均为公式推导使用的中间变量，且中间变量
Figure GDA00022228003300000918
和
Figure GDA00022228003300000919
的计算函数式如式(4)所示；
Figure GDA0002222800330000093
式(4)中，A、B为公式推导使用的中间变量，
Figure GDA00022228003300000920
表示t0时刻平均轨道要素偏差到tf时刻平均轨道要素偏差的一阶传递矩阵，
Figure GDA00022228003300000921
表示t0时刻轨道要素偏差的无量纲化矩阵，
Figure GDA00022228003300000922
表示tf时刻轨道要素偏差的无量纲化矩阵，
Figure GDA00022228003300000923
表示t0时刻平均轨道要素偏差到tf时刻平均轨道要素偏差的二阶传递张量，Q(θf)表示tf时刻无量纲轨道要素偏差到相对状态的二阶转换张量，T(θf)表示tf时刻以纬度幅角为自变量的无量纲化坐标到以时间为自变量的量纲化坐标的转换矩阵，T-1(θf)表示T(θf)的逆矩阵，Π表示相对状态分量次序变换矩阵，Σ(θf)表示tf时刻轨道要素偏差到相对状态的转换矩阵，
Figure GDA00022228003300000924
表示平均轨道要素偏差到吻切轨道要素偏差的转换矩阵，
Figure GDA00022228003300000925
为Kronecker张量积。本实施例中，
Figure GDA00022228003300000926
等带上标“bar”的矩阵或张量表示与平均轨道根数对应的值，即将主航天器对应时刻的平均轨道根数带入矩阵D，Γ，H表达式计算获得的值，反之没有上标“bar”的矩阵或张量表示与吻切轨道根数对应的值。式(4)中有对应初始时刻t0及任意时刻tf的量，分别带入t0时刻或tf时刻下主航天器的轨道参数，即可计算对应的矩阵及张量。将θf＝θ0带入式(4)，即可计算
Figure GDA0002222800330000102
及
Figure GDA0002222800330000103
其中
Figure GDA0002222800330000104
In为n维单位矩阵。需要说明的是，前述矩阵T，T-1，Π，Γ，P及张量Q、H均为已知矩阵，其详细表达式参见文献:SenguptaP,Vadali S R,Alfriend K T.Second-order state transition for relative motionnear perturbed,elliptic orbits[J].Celestial Mechanics and DynamicalAstronomy,2006,97(2):101-129。前述矩阵Σ，
Figure GDA0002222800330000105
D同样也均为已知矩阵，其详细表达式参见文献:Gim D W,Alfriend K T.State Transition Matrix of Relative Motion forthe Perturbed Noncircular Reference Orbit[J].Journal of Guidance,Control,andDynamics,2003,26(6):956–971。
本实施例中，步骤3)计算用于解析预报卫星相对运动状态偏差的一阶及二阶状态转移张量的函数表达式如式(5)所示；
Figure GDA0002222800330000101
式(5)中，dxf为终端时刻的相对运动状态偏差，
Figure GDA0002222800330000106
为用于解析预报卫星相对运动状态偏差的一阶状态转移张量，dx0为初始的相对运动状态偏差，
Figure GDA0002222800330000107
为用于解析预报卫星相对运动状态偏差的二阶状态转移张量，
Figure GDA0002222800330000108
为
Figure GDA0002222800330000109
的第i行第j列元素，
Figure GDA00022228003300001010
为
Figure GDA00022228003300001011
的第i行第j列元素，
Figure GDA00022228003300001012
为
Figure GDA00022228003300001013
的第i行k列j维元素，
Figure GDA00022228003300001014
是从t0时刻到tf时刻的状态转移矩阵，Φ(tf,t0)是从t0时刻到tf时刻的状态转移矩阵，Ψ(tf,t0)为从t0时刻到tf时刻的二阶状态转移张量，
Figure GDA00022228003300001015
为初始时刻t0从卫星相对主卫星的标称相对运动状态δx0的第k个分量，
Figure GDA00022228003300001016
为Kronecker张量积。
本实施例中，步骤5)中获得的卫星相对运动状态偏差的均值及协方差矩阵的解析非线性传播方程的函数表达式如式(6)所示；
Figure GDA0002222800330000111
式(6)中，mi(tf)为卫星相对状态偏差的均值矩阵m(tf)的第i行元素，mj(tf)为卫星相对状态偏差的均值矩阵m(tf)的第j行元素，Cij(tf)为卫星相对状态偏差的矩阵，
Figure GDA0002222800330000113
为用于解析预报卫星相对运动状态偏差的一阶状态转移张量
Figure GDA0002222800330000114
的第i行第a列元素，
Figure GDA0002222800330000115
为用于解析预报卫星相对运动状态偏差的一阶状态转移张量
Figure GDA0002222800330000116
的第j行第b列元素，
Figure GDA0002222800330000117
为从t0时刻到tf时刻的二阶状态转移张量
Figure GDA0002222800330000118
的第i行a列b维元素，
Figure GDA0002222800330000119
为从t0时刻到tf时刻的二阶状态转移张量
Figure GDA00022228003300001110
的第j行b列c维元素，
Figure GDA00022228003300001111
为从t0时刻到tf时刻的二阶状态转移张量
Figure GDA00022228003300001112
的第i行b列c维元素，
Figure GDA00022228003300001113
为从t0时刻到tf时刻的二阶状态转移张量
Figure GDA00022228003300001114
的第j行c列d维元素，
Figure GDA00022228003300001115
分别为初始相对状态偏差dx0的前四阶矩。
本实施例中，初始相对状态偏差dx0的前四阶矩的计算函数表达式如式(7)所示；
Figure GDA0002222800330000112
式(7)中，
Figure GDA00022228003300001116
分别为初始相对状态偏差dx0的前四阶矩，
Figure GDA00022228003300001117
分别为初始相对状态偏差的均值矩阵m0的a,b,c,d分量元素，
Figure GDA00022228003300001118
分别为初始相对状态偏差的协方差矩阵C0的ab分量元素、bc分量元素、cd分量元素、ad分量元素。本实施例中，上标a,b,c,d为矩阵及张量对应的分量元素。对高斯分布，概率密度函数仅有均值及协方差矩阵确定，因此终端相对状态偏差的概率密度函数可表示为pg(dxf；mf,Cf)。需要说明的是，即使初始相对状态偏差为高斯分布，经非线性传播后，终端偏差一般也为非高斯的。因为仅用均值和协方差矩阵无法完整描述非高斯偏差的概率密度函数，所以本发明用高斯和模型来预报非高斯分布偏差的概率密度函数。
本实施例中，步骤6)的详细步骤包括：
6.1)用N个高斯分布概率密度函数去逼近初始相对状态偏差的概率密度函数p(t0,dx0)获得每一个子高斯分布的权重、均值和协方差矩阵；
6.2)采用卫星相对运动状态偏差的均值及协方差矩阵的解析非线性传播方程，将每个子高斯分布的均值mi(t0)及协方差矩阵Ci(t0)预报到终端时刻tf，预报过程中子高斯分布的权重保持不变，获得终端时刻每个子高斯的均值mi(tf)及协方差矩阵Ci(tf)，其中i＝1,2,…,N，N为高斯分布概率密度函数的数量；
6.3)针对每一个第i个子高斯分布，其中i＝1,2,…,N，N为高斯分布概率密度函数的数量，根据其权重系数ωi、均值mi(tf)及协方差矩阵Ci(tf)，根据式(8)和式(9)所示获得卫星非线性相对运动状态偏差的概率密度函数的解析传播方程计算卫星相对状态偏差的均值及协方差矩阵以及卫星偏差的概率密度函数p(tf,dxf)；
Figure GDA0002222800330000121
Figure GDA0002222800330000122
式(8)和(9)中，m(tf),C(tf)为卫星相对状态偏差dxf的均值及协方差矩阵，ωi表示第i个子高斯分布的权重系数，N为高斯分布概率密度函数的数量，ω为所有高斯分布的权重系数之和，mi(tf)为tf时刻第i个子高斯分布的均值，Ci(tf)为tf时刻第i个子高斯分布的协方差矩阵，
Figure GDA0002222800330000124
为mi(tf)的转置矩阵，mT(tf)为m(tf)的转置矩阵，p(tf,dxf)为卫星偏差的概率密度函数，pg(dxf；mi(tf),Ci(tf))表示tf时刻第i个子高斯分布的概率密度函数。
本实施例中，步骤6.1)中逼近初始相对状态偏差的概率密度函数p(t0,dx0)采用的逼近公式如式(10)所示；
Figure GDA0002222800330000123
式(10)中，p(t0,dx0)为初始相对状态偏差的概率密度函数，dx0为初始相对状态偏差，pg(dx0；mi,Ci)表示第i个子高斯分布的概率密度函数，mi和Ci分别为相应的均值和协方差矩阵，ωi表示第i个子高斯分布的权重系数。
本实施例中，第i个子高斯分布的权重系数ωi的计算函数表达式如式(11)所示；
Figure GDA0002222800330000131
式(11)中，ωi表示第i个子高斯分布的权重系数，m0、C0为初始相对状态偏差的均值及协方差矩阵，
Figure GDA0002222800330000132
及
Figure GDA0002222800330000133
为将一维标准正态分布pg(x；0,1)等方差地分割为N个子正态分布的权重、均值及标准差，diag{a}表示以a向量为对角线元素的矩阵，V和Λ分别为对C0进行特征值分解所得的特征向量矩阵与特征值矩阵，λ1,λ2,...,λn为C0的特征值，λj与υj分别为沿选定的j方向的特征值与特征向量，即Λ的第j行j列元素和V的第j列向量，j选择为初始偏差协方差矩阵中径向位置或横向速度的对应值。需要说明的是，
Figure GDA0002222800330000134
及
Figure GDA0002222800330000135
可通过优化方法求解，也可直接调用现有的数据库获得，见文献Vittaldev,V.,and Russell,R.P.,Space Object Collision Probability Using Multidirectional GaussianMixture Models[J].Journal of Guidance,Control,and Dynamics,2016,39(9):2163-2169。
综上所述，本实施例解析的卫星非线性相对运动偏差传播分析方法为一种考虑J2摄动、解析的卫星非线性相对运动偏差传播分析方法，利用编队卫星的初始相对运动状态均值及协方差矩阵，计算相对状态及相对状态偏差的状态转移矩阵及张量，预报相对运动状态偏差的均值、协方差矩阵、及概率密度函数。其中，可根据初始偏差及终端偏差是否为高斯分布，选择是否采用高斯和模型。该方法考虑了J2摄动项和二阶非线性项、可用于相距较远航天器的长时间、高精度相对运动状态偏差传播，具有设计方法正确合理，方法解析、计算速度快，对实际工程任务的适用性好等优点。
参见图2所示终端误差分布对比可知，和精确的蒙特卡洛仿真相比，本方法用的高斯和模型的所有子高斯分布很好的包围了的蒙特卡洛仿真的样本点，说明本方法精确可靠。进一步由图3可知，本方法用的高斯和模型对均值和协方差的预报与蒙特卡洛仿真一致，其概率密度等高线与蒙特卡洛仿真的样本点分布一致，且很好的表征了终端偏差的非高斯分布特性，说明本方法在非线性、非高斯偏差均值、协方差矩阵及概率密度函数的传播分析中具有较高的精度。
实施例二：
本实施例与实施例一基本相同，其主要不同点为：本实施例中，步骤3)中初始及终端偏差均为高斯分布，因此不再进行步骤6)的计算，仅需要采用步骤5)对非线性相对状态偏差的均值及协方差矩阵进行解析预报，终端偏差的概率密度函数由其均值及协方差矩阵即可确定。本实施例仅需要预报相对状态偏差的均值及协方差矩阵，不需要采用高斯和模型预报概率密度函数，方法解析、计算效率高，实施简单；计算效率高。参见图4所示终端误差均值及协方差矩阵的预报结果对比可知，本本实施例方法用的高斯和模型对均值和协方差矩阵的预报与蒙特卡洛仿真一致，相比已有的线性方法，本方法精度显著提高。
实施例三：
本实施例与实施例一基本相同，其主要不同点为：本实施例中，不再进行步骤2)的计算，将步骤3)中(6)式的
Figure GDA0002222800330000141
矩阵直接用步骤1)中(1)式的
Figure GDA0002222800330000142
矩阵替换。因为卫星绝对状态偏差可作为一个相对状态，采用(1)式进行预报，所以本实施例可用于解析地预报卫星绝对状态偏差的均值、协方差矩阵及概率密度函数。本实施例考虑了J2摄动项和二阶非线性项、可用于单个航天器长时间、高精度的绝对状态偏差传播，方法解析、计算效率高。
以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (9)

1.一种解析的卫星非线性相对运动偏差传播分析方法，其特征在于实施步骤包括：
1)为编队飞行的两卫星指定主卫星与从卫星，并输入初始时刻的参考卫星绝对轨道状态、两卫星标称相对运动状态、输入初始相对运动状态偏差的概率密度函数；
2)基于初始时刻的参考卫星绝对轨道状态、两卫星标称相对运动状态，根据考虑J2摄动的非线性相对运动方程计算用于解析预报卫星相对运动状态的一阶及二阶状态转移张量；
3)基于卫星相对运动状态的一阶及二阶状态转移张量，计算用于解析预报卫星相对运动状态偏差的一阶及二阶状态转移张量；
4)判断初始及终端时刻的卫星偏差是否均为高斯分布，如果均为高斯分布则跳转执行步骤5)，否则，跳转执行步骤6)；
5)采用协方差分析方法，获得卫星相对运动状态偏差的均值及协方差矩阵的解析非线性传播方程，跳转执行步骤7)；
6)采用高斯和模型，按步骤6.1)～6.3)获得卫星非线性相对运动状态偏差的概率密度函数的解析传播方程，然后跳转执行步骤7)；
6.1)用N个子高斯分布概率密度函数去逼近初始相对状态偏差的概率密度函数p(t0,dx0)获得每一个子高斯分布的权重ωi、均值mi(t0)和协方差矩阵Ci(t0)；
6.2)采用卫星相对运动状态偏差的均值及协方差矩阵的解析非线性传播方程，将每个子高斯分布的均值mi(t0)及协方差矩阵Ci(t0)预报到终端时刻tf，预报过程中子高斯分布的权重ωi保持不变，获得卫星终端时刻每个子高斯的均值mi(tf)及协方差矩阵Ci(tf)，其中i＝1,2,…,N，N为子高斯分布概率密度函数的数量；
6.3)针对第i个子高斯分布，其中i＝1,2,…,N，N为子高斯分布概率密度函数的数量，根据其权重系数ωi、均值mi(tf)及协方差矩阵Ci(tf)，获得卫星非线性相对运动状态偏差的解析传播均值为
Figure FDA0002289590700000011
协方差矩阵为
Figure FDA0002289590700000012
概率密度函数为
Figure FDA0002289590700000013
其中
Figure FDA0002289590700000014
dxf为卫星终端相对运动状态偏差，pg(dxf；mi(tf),Ci(tf))表示tf时刻第i个子高斯分布的概率密度函数；
7)输出终端时刻两卫星相对运动状态偏差的均值、协方差矩阵及概率密度函数。
2.根据权利要求1所述的解析的卫星非线性相对运动偏差传播分析方法，其特征在于，步骤1)中输入初始时刻的参考卫星绝对轨道状态具体是指输入主卫星初始时刻t0的绝对轨道状态E0＝[a0,e0,i0,Ω0,ω0,f0]，其中，E0为主卫星初始时刻t0的绝对轨道状态，a0为主卫星初始时刻t0的半长轴，e0为主卫星初始时刻t0的偏心率，i0为主卫星初始时刻t0的轨道倾角，Ω0为主卫星初始时刻t0的升交点赤经，ω0为主卫星初始时刻t0的近拱点角距，f0为主卫星初始时刻t0的真近点角；且输入初始时刻t0从卫星相对主卫星的标称相对运动状态为δx(t0)，δx(t0)在主卫星当地轨道坐标系中表示，该坐标系原点为主卫星质心，x轴沿主卫星地心矢径方向，z轴沿轨道面法向，y轴与x、z轴构成右手坐标系；输入初始相对运动状态偏差为dx(t0)及其概率密度函数为p(t0,dx0)，若初始偏差为高斯分布，则其概率密度函数为pg(dx0；m0,C0)，m0为初始相对状态偏差的均值矩阵，C0为初始相对状态偏差的协方差矩阵。
3.根据权利要求2所述的解析的卫星非线性相对运动偏差传播分析方法，其特征在于，步骤2)中考虑J2摄动的非线性相对运动方程的函数表达式如式(1)所示；
Figure FDA0002289590700000021
式(1)中，δx(tf)表示时刻tf从卫星相对主卫星的标称相对运动状态，Φ(tf,t0)是从t0时刻到tf时刻的状态转移矩阵，δx(t0)表示初始时刻t0从卫星相对主卫星的标称相对运动状态，Ψ(tf,t0)为从t0时刻到tf时刻的二阶状态转移张量，T(tf)为以纬度幅角θf为自变量的无量纲化坐标到以时间tf为自变量的量纲化坐标的转换矩阵，
Figure FDA0002289590700000022
为从θ0纬度幅角到θf纬度幅角的一阶状态转移张量，T(t0)为以纬度幅角θ0为自变量的无量纲化坐标到以时间t0为自变量的量纲化坐标的转换矩阵，
Figure FDA0002289590700000023
为从θ0纬度幅角到θf纬度幅角的二阶状态转移张量，T-1(t0)为转换矩阵T(t0)的逆矩阵，
Figure FDA0002289590700000024
为Kronecker张量积。
4.根据权利要求3所述的解析的卫星非线性相对运动偏差传播分析方法，其特征在于，步骤2)中计算用于解析预报卫星相对运动状态的一阶状态转移张量如式(2)所示，计算用于解析预报卫星相对运动状态的二阶状态转移张量如式(3)所示；
Figure FDA0002289590700000025
Figure FDA0002289590700000026
式(2)和式(3)中，
Figure FDA0002289590700000027
为用于解析预报卫星相对运动状态的从θ0纬度幅角到θf纬度幅角的一阶状态转移张量，
Figure FDA0002289590700000031
表示θf＝θ0时的中间变量
Figure FDA0002289590700000032
表示
Figure FDA0002289590700000033
的逆矩阵，
Figure FDA0002289590700000034
表示θf＝θ0时的中间变量
Figure FDA0002289590700000035
为用于解析预报卫星相对运动状态的从θ0纬度幅角到θf纬度幅角的二阶状态转移张量，
Figure FDA0002289590700000036
为Kronecker张量积，
Figure FDA0002289590700000037
和
Figure FDA0002289590700000038
均为公式推导使用的中间变量，且中间变量
Figure FDA0002289590700000039
和
Figure FDA00022895907000000310
的计算函数式如式(4)所示；
Figure FDA00022895907000000311
式(4)中，A、B为公式推导使用的中间变量，
Figure FDA00022895907000000312
表示t0时刻平均轨道要素偏差到tf时刻平均轨道要素偏差的一阶传递矩阵，
Figure FDA00022895907000000313
表示t0时刻轨道要素偏差的无量纲化矩阵，
Figure FDA00022895907000000314
表示tf时刻轨道要素偏差的无量纲化矩阵，
Figure FDA00022895907000000315
表示t0时刻平均轨道要素偏差到tf时刻平均轨道要素偏差的二阶传递张量，Q(θf)表示tf时刻无量纲轨道要素偏差到相对状态的二阶转换张量，T(θf)表示tf时刻以纬度幅角为自变量的无量纲化坐标到以时间为自变量的量纲化坐标的转换矩阵，T-1(θf)表示T(θf)的逆矩阵，Π表示相对状态分量次序变换矩阵，Σ(θf)表示tf时刻轨道要素偏差到相对状态的转换矩阵，
Figure FDA00022895907000000316
表示平均轨道要素偏差到吻切轨道要素偏差的转换矩阵，
Figure FDA00022895907000000317
为Kronecker张量积。
5.根据权利要求4所述的解析的卫星非线性相对运动偏差传播分析方法，其特征在于，步骤3)计算用于解析预报卫星相对运动状态偏差的一阶及二阶状态转移张量的函数表达式如式(5)所示；
Figure FDA00022895907000000318
式(5)中，dxf为终端时刻的相对运动状态偏差，
Figure FDA00022895907000000319
为用于解析预报卫星相对运动状态偏差的一阶状态转移张量，dx0为初始的相对运动状态偏差，
Figure FDA00022895907000000320
为用于解析预报卫星相对运动状态偏差的二阶状态转移张量，
Figure FDA00022895907000000321
为
Figure FDA00022895907000000322
的第i行第j列元素，为
Figure FDA00022895907000000324
的第i行第j列元素，
Figure FDA00022895907000000325
为
Figure FDA00022895907000000326
的第i行k列j维元素，
Figure FDA00022895907000000327
是从t0时刻到tf时刻的状态转移矩阵，
Figure FDA00022895907000000328
为从t0时刻到tf时刻的二阶状态转移张量，
Figure FDA00022895907000000329
为初始时刻t0从卫星相对主卫星的标称相对运动状态δx0的第k个分量，
Figure FDA0002289590700000041
为Kronecker张量积。
6.根据权利要求1所述的解析的卫星非线性相对运动偏差传播分析方法，其特征在于，步骤5)中获得的卫星相对运动状态偏差的均值及协方差矩阵的解析非线性传播方程的函数表达式如式(6)所示；
Figure FDA0002289590700000042
式(6)中，mi(tf)为卫星相对状态偏差的均值矩阵m(tf)的第i行元素，mj(tf)为卫星相对状态偏差的均值矩阵m(tf)的第j行元素，Cij(tf)为卫星相对状态偏差的矩阵，
Figure FDA0002289590700000043
为用于解析预报卫星相对运动状态偏差的一阶状态转移张量
Figure FDA0002289590700000044
的第i行第a列元素，
Figure FDA0002289590700000045
为用于解析预报卫星相对运动状态偏差的一阶状态转移张量
Figure FDA0002289590700000046
的第j行第b列元素，
Figure FDA0002289590700000047
为从t0时刻到tf时刻的二阶状态转移张量
Figure FDA0002289590700000048
的第i行a列b维元素，
Figure FDA0002289590700000049
为从t0时刻到tf时刻的二阶状态转移张量
Figure FDA00022895907000000410
的第j行b列c维元素，
Figure FDA00022895907000000411
为从t0时刻到tf时刻的二阶状态转移张量
Figure FDA00022895907000000412
的第i行b列c维元素，
Figure FDA00022895907000000413
为从t0时刻到tf时刻的二阶状态转移张量
Figure FDA00022895907000000414
的第j行c列d维元素，
Figure FDA00022895907000000415
分别为初始相对状态偏差dx0的前四阶矩。
7.根据权利要求6所述的解析的卫星非线性相对运动偏差传播分析方法，其特征在于，初始相对状态偏差dx0的前四阶矩的计算函数表达式如式(7)所示；
Figure FDA00022895907000000416
式(7)中，
Figure FDA00022895907000000417
分别为初始相对状态偏差dx0的前四阶矩，
Figure FDA00022895907000000418
分别为初始相对状态偏差的均值矩阵m0的a,b,c,d分量元素，
Figure FDA00022895907000000419
分别为初始相对状态偏差的协方差矩阵C0的ab分量元素、bc分量元素、cd分量元素、ad分量元素。
8.根据权利要求1所述的解析的卫星非线性相对运动偏差传播分析方法，其特征在于，步骤6.1)中逼近初始相对状态偏差的概率密度函数p(t0,dx0)采用的逼近公式如式(10)所示；
Figure FDA0002289590700000051
式(10)中，p(t0,dx0)为初始相对状态偏差的概率密度函数，dx0为初始相对状态偏差，pg(dx0；mi,Ci)表示第i个子高斯分布的概率密度函数，mi和Ci分别为相应的均值和协方差矩阵，ωi表示第i个子高斯分布的权重系数。
9.根据权利要求8所述的解析的卫星非线性相对运动偏差传播分析方法，其特征在于，第i个子高斯分布的权重系数ωi的计算函数表达式如式(11)所示；
Figure FDA0002289590700000052
式(11)中，ωi表示第i个子高斯分布的权重系数，m0、C0为初始相对状态偏差的均值及协方差矩阵，
Figure FDA0002289590700000053
及
Figure FDA0002289590700000054
为将一维标准正态分布pg(x；0,1)等方差地分割为N个子正态分布的权重、均值及标准差，diag{a}表示以a向量为对角线元素的矩阵，V和Λ分别为对C0进行特征值分解所得的特征向量矩阵与特征值矩阵，λ1,λ2,...,λn为C0的特征值，λj与υj分别为沿选定的j方向的特征值与特征向量，即Λ的第j行j列元素和V的第j列向量，j选择为初始偏差协方差矩阵中径向位置或横向速度的对应值。
CN201710287256.9A 2017-04-27 2017-04-27 一种解析的卫星非线性相对运动偏差传播分析方法 Active CN106970643B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710287256.9A CN106970643B (zh) 2017-04-27 2017-04-27 一种解析的卫星非线性相对运动偏差传播分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710287256.9A CN106970643B (zh) 2017-04-27 2017-04-27 一种解析的卫星非线性相对运动偏差传播分析方法

Publications (2)

Publication Number Publication Date
CN106970643A CN106970643A (zh) 2017-07-21
CN106970643B true CN106970643B (zh) 2020-03-27

Family

ID=59334055

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710287256.9A Active CN106970643B (zh) 2017-04-27 2017-04-27 一种解析的卫星非线性相对运动偏差传播分析方法

Country Status (1)

Country Link
CN (1) CN106970643B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107741240B (zh) * 2017-10-11 2020-11-24 成都国卫通信技术有限公司 一种适用于动中通的组合惯导系统自适应初始对准方法
CN108490966B (zh) * 2018-01-31 2021-02-05 中国人民解放军国防科技大学 基于微分代数的静止轨道摄动相对轨迹高阶制导方法
CN109255096B (zh) * 2018-07-25 2022-10-04 西北工业大学 一种基于微分代数的地球同步卫星轨道不确定演化方法
CN110053788B (zh) * 2019-03-15 2022-05-13 中国西安卫星测控中心 一种考虑复杂摄动的星座长期保持控制频次估计方法
CN110442831B (zh) * 2019-07-31 2023-03-24 中国人民解放军国防科技大学 基于非线性偏差演化的空间非合作目标天基搜索方法
CN112800651B (zh) * 2021-01-27 2022-10-18 重庆大学 一种等体积超球体筛点法的结构可靠度分析方法
CN113836697B (zh) * 2021-08-26 2024-01-23 北京理工大学 基于动态混合高斯的空间引力波探测构型稳定性演化方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6708116B2 (en) * 2001-11-13 2004-03-16 Analytical Graphics, Inc. Method and apparatus for orbit determination
CN105069311B (zh) * 2015-08-24 2018-06-12 哈尔滨工业大学 一种远程火箭发射初态误差传播估计方法
CN106297422B (zh) * 2016-10-09 2019-01-08 中国人民解放军国防科学技术大学 一种基于最小二乘的非线性相对轨迹预报方法

Also Published As

Publication number Publication date
CN106970643A (zh) 2017-07-21

Similar Documents

Publication Publication Date Title
CN106970643B (zh) 一种解析的卫星非线性相对运动偏差传播分析方法
Luo et al. A review of uncertainty propagation in orbital mechanics
Shen et al. Novel variable structure measurement system with intelligent components for flight vehicles
Li et al. A machine learning-based approach for improved orbit predictions of LEO space debris with sparse tracking data from a single station
Yang et al. Comparison of unscented and extended Kalman filters with application in vehicle navigation
Sabet et al. Optimal design of the own ship maneuver in the bearing-only target motion analysis problem using a heuristically supervised extended Kalman filter
Jones et al. A labeled multi-Bernoulli filter for space object tracking
CN109508023B (zh) 俯角参考跟踪系统
Wang et al. Onboard satellite visibility prediction using metamodeling based framework
Zeng et al. Potential hop reachable domain over surfaces of small bodies
Sanwale et al. Aerodynamic Parameters Estimation Using Radial Basis Function Neural Partial Differentiation Method.
Lou et al. Consider unobservable uncertain parameters using radio beacon navigation during Mars entry
Chai et al. An interactive fuzzy physical programming for solving multiobjective skip entry problem
Jingsen et al. Integrating extreme learning machine with Kalman filter to bridge GPS outages
Ma et al. Jet transport particle filter for attitude estimation of tumbling space objects
Paul et al. Advanced ensemble modeling method for space object state prediction accounting for uncertainty in atmospheric density
Buzhardt et al. Controlled density transport using Perron Frobenius generators
Zhou et al. Inverse simulation system for manual-controlled rendezvous and docking based on artificial neural network
Sabatini et al. Carrier-phase GNSS Attitude Determination and Control for Small Unmanned Aerial Vehicle Applications
Huang et al. Uncertainty analysis of reachable set for planetary entry using polynomial chaos
Ridderhof et al. Planetary entry in a randomly perturbed atmosphere
Blanke et al. Structural analysis—a case study of the Rømer satellite
Aguiar et al. Kalman filtering for differential drive robots tracking
de la Mota et al. Data-driven probabilistic methodology for aircraft conflict detection under wind uncertainty
Qi et al. Hybrid Smith predictor and phase lead based divergence compensation for hardware-in-the-loop contact simulation with measurement delay

Legal Events

Date Code Title Description
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