CN110779531A - 一种天基仅测角差分进化一次精确定轨方法 - Google Patents

一种天基仅测角差分进化一次精确定轨方法 Download PDF

Info

Publication number
CN110779531A
CN110779531A CN201910870972.9A CN201910870972A CN110779531A CN 110779531 A CN110779531 A CN 110779531A CN 201910870972 A CN201910870972 A CN 201910870972A CN 110779531 A CN110779531 A CN 110779531A
Authority
CN
China
Prior art keywords
vector
observation
orbit
target
angle
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.)
Granted
Application number
CN201910870972.9A
Other languages
English (en)
Other versions
CN110779531B (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.)
Beijing University of Aeronautics and Astronautics
Original Assignee
Beijing 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 Beijing University of Aeronautics and Astronautics filed Critical Beijing University of Aeronautics and Astronautics
Priority to CN201910870972.9A priority Critical patent/CN110779531B/zh
Publication of CN110779531A publication Critical patent/CN110779531A/zh
Application granted granted Critical
Publication of CN110779531B publication Critical patent/CN110779531B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Astronomy & Astrophysics (AREA)
  • Automation & Control Theory (AREA)
  • General Physics & Mathematics (AREA)
  • Navigation (AREA)

Abstract

本发明针对太空非合作目标精定轨问题，提出了一种天基仅测角差分进化一次精确定轨方法。该方法基于最小二乘批处理定轨思路，使用差分进化算法直接优化定轨目标函数，不仅避免了直接最小二乘批处理方法求解过程可能出现的不收敛或收敛错误情形，而且，不需要一般精定轨要求的二次定轨过程，可以达到一次精确定轨目的。该方法基于轨道动力学方程描述形式构造目标轨道参数向量，在仅使用角度观测数据的情况下，精确求取目标在特定时刻的轨道空间位置和速度。仿真实验表明，使用精度为5角秒的测角序列，有效观测距离在5000公里范围内，该方法的定轨精度可达百米量级，可用于卫星定轨任务。

Description

一种天基仅测角差分进化一次精确定轨方法
技术领域
本发明属于卫星定轨技术领域，尤其涉及一种天基差分进化的仅测角一次精确定轨方法。
背景技术
太空存在大量空间碎片、故障航天器等非合作目标，这些目标对航天器在太空中的飞行造成了很大的安全隐患。因此，对非合作目标的监视在空间任务中有着重要意义。空间监视可分为地基监视和天基监视两类。相较于地基监(GBSS)，天基监视不受观测站位置和数量约束，可以根据任务需求获取更多观测数据。天基监视中，光学成像是针对非合作目标的主要观测方式。通过光学成像方式观测目标，可获取目标相对于观测航天器之间的空间角度信息，继而由连续观测目标空间角度的变化得到非合作目标特定时刻的轨道参数。
仅测角天基目标定轨方法使用目标的空间角度观测序列和天基观测航天器轨道信息估计目标轨道信息，其求解基于观测方程和卫星轨道动力学方程约束，使用观测角度序列估计目标在特定时刻的轨道参数，属于非线性状态估计问题。定轨的最终目标是得到非合作目标的精确轨道参数。目前已有的精确定轨方法分两步完成：初定轨和轨道优化，其中，初定轨的结果距离真值有较大偏差，仅用于为轨道优化提供合适的初值；轨道优化在初值基础上通过迭代求解得到精确轨道参数。
初定轨经典方法主要基于Laplace、Gauss和双r迭代法得到目标轨道近似解，这些方法基于位置矢量的空间几何关系，使用测角信息，建立方程求解轨道的位置和速度信息。此外，进化方法也被用于初定轨求解当中，在轨道可能的状态空间中搜寻使推测角度序列最接近真实测量数据的轨道初值，相较于Gauss、Laplace传统初定轨，进化方法避免了传统解算中矩阵奇异、迭代不收敛的问题。Ansalone等首先使用遗传算法搜索Lambert问题解空间，继而通过计算轨道运动学参数，实现了目标初定轨；Hideaki Hinagawa等同样基于遗传算法求解Lambert问题，针对实际地球同步轨道数据进行了初定轨的有效性验证。Li等使用遗传算法选取了和L.Ansalone不同的搜索变量和轨道参数计算策略，减少了观测数据噪声对定轨精度的影响，同年，Li等在遗传算法初定轨基础上，基于抗差估计理论选取了LTS(LeastTrimmed Square)方法设计评价函数，并通过与其他评价函数比较，验证了选取的评价函数可以增加定轨对观测误差和出格点(outlier)的容忍度；除遗传算法外，Li等结合分布估计算法(EDA)和差分进化算法(DE)算法，通过寻找最优解最大概率分布的点确定轨道参数，减少了观测误差对定轨结果的影响。
通常精定轨需要两个计算过程，一是，初定轨确定轨道近似参数；二是，轨道优化，在初定轨基础上进行，以初定轨结果为初值进行迭代计算，得到目标轨道的精确数据。目前，轨道优化方法通常使用最小二乘批处理算法。最小二乘批处理方法基于观测角度序列残差平方和最小的原则，采用迭代求解的方式得到目标航天器的精确轨道参数。已有的轨道优化方法均对初值选取有较高的要求，偏离真实轨道数据较远的初值会造成定轨结果不准确甚至发散。此外，现有轨道优化方法在运算过程中均对轨道转移矩阵进行了线性化近似，引入了计算误差。轨道优化方法在求解过程中还涉及了正规矩阵求逆问题，这同样会引入计算误差，在正规矩阵条件数较大时，计算结果会严重偏离真实值，甚至无法收敛。
发明内容
为解决上述问题，本文基于最小二乘批处理轨道优化思路，设计了差分进化仅测角一次精定轨方法。差分进化定轨相较现有的轨道优化方法在目标状态解算环节使用了不同策略。差分进化方法直接在基于轨道动力学方程构造的参数空间中搜索最优解，估计目标轨道参数残差平方和的全局最小状态，省略了初定轨需求，避免了运算中的线性化误差，并且不需要对正规矩阵进行计算，避免了计算过程引入的误差，可以提高定轨精度。
另外，相较于已有进化算法在定轨中的应用，本发明考虑了轨道摄动影响，基于轨道动力学方程描述方式，将目标位置和速度设为搜索变量，并在最小二乘法思路上对评价函数进行了进一步设计，提高了观测误差干扰下的定轨精度。本文使用STK生成的轨道数据进行仿真验证，结果证明了算法的有效性。
采用下述技术方案：一种用于太空非合作目标的天基差分进化仅测角一次精定轨方法，其特征在于，包括如下步骤：
S1：将卫星差t0时刻的位置速度状态作为参数，记为
Figure BDA0002202810440000021
并设置差分化进化算法参数，所述参数包括变量维数D、种群个数NP、变异缩放因子F；变异率Cr；最大迭代代数gmax，其中，D为6；
S2：根据步骤S1中设置的变量维数D、种群个数NP进行种群初始化，生成初始种群，种群个体的目标星轨道参数预测值为Xk,g，k＝1…NP；
S3：对步骤S2生成的初始种群，计算以种群个体Xk,g为初始状态时，在各观测时刻观测方向与估计方向夹角θi(i＝1,…,n，n为观测点个数)，将夹角θi作为i时刻的估计误差，将各时刻误差序列表示为向量形式：r＝[θ0,θ1,…,θn]T。
S4：通过最小二乘法得到估计值残差定义，设置估计值残差为RSM，通过式(1)对步骤S2生成种群的每个个体进行残差值计算；
RMS＝rTr (1)
其中r＝[θ0,θ1,…,θn]T。
S5：检查并剔除步骤S4中误差向量r中的异常值，若误差向量r的第i项|ri|≥3σ，其中，即看作异常值；
其中，σ2为r各分量平方的均值，
S6：使用步骤S4得到的残差值RMS和步骤S5得到的异常值计算评价函数，得到Xk,g的评价函数值ffitness(X)：评价函数如式(2)所示，
Figure BDA0002202810440000032
其中，m为异常值个数，n-m为非异常值个数，routlier为差值向量r中m个异常值组成的向量；
S7：变异，通过比较步骤S6中获得的Xk,g的评价函数值ffitness(X)，得出评价函数最优状态向量Xbest,g，并通过最优状态向量Xbest,g进行变异得到交叉变异状态向量vk,g；
S8：交叉，通过交叉变异状态向量vk,g和原状态向量Xk,g中不同位置元素，产生新的候选状态向量uk,g；
S9：选择，通过对比候选状态向量uk,g和原状态向量Xk,g的评价函数，生成下一代种群的状态向量Xk,g+1；
S10：终止条件，设第g代最佳状态向量Xbest,g对应的评价函数为fbest,g，设最大进化代数为gmax；当达到以下条件之一时，差分进化计算过程终止：
(1)fbest,g-fbest,g-1＜δ，δ正整数
(2)g＞gmax
取终止计算时的Xbest,g作为定轨结果。
优选地，所述步骤S2中种群的初始化过程具体为：
设有效观测距离为d，观测平台轨道参数为
Figure BDA0002202810440000042
在状态空间中，rmin表示航天器与地心距离的最小值，根据rmin和目标偏心率范围可以计算出速度最大值vmax，即坐标轴各方向的速度取值范围为[-vmax,vmax]；
根据公式(3)进行种群初始化过程，
Figure BDA0002202810440000041
其中，状态向量上界bU＝(xT+d,yT+d,zT+d,vmax,vmax,vmax)T，状态变量下界bL＝(xT-d,yT-d,zT-d,-vmax,-vmax,-vmax)T，Xj,k,g为第g代第k个个体的第j项；P为种群，表示所有的状态向量，g为当前迭代代数；D为变量维数。
优选地，步骤S3中的观测方向与估计方向夹角θi通过公式(4)确定：
Figure BDA0002202810440000051
其中，设ti时刻观测平台坐标系下观测矢量为Li，估计矢量为
Figure BDA0002202810440000052
优选地，观测矢量L根据公式(5)：
L＝ρ[cosαcosβ,sinαcosβ,sinβ]T,ρ＞0 (5)
其中，ρ表示观测平台与目标星间的距离，因距离大小不影响矢量方向，计算中，ρ可取任意正整数；α、β为实际观测的俯仰角和方位角。
优选地，α、β根据公式(6)确定，其推算过程如下：
当t0时刻估计初值为X0时，可以得到ti时刻目标星在观测平台坐标系下的位置rs＝[xs,ys,zs]T时，对应的观测角的可通过式(6)计算：
Figure BDA0002202810440000053
其中，xs,ys和zs为目标在观测平台坐标系下的空间位置，其中，H＝[αβ]T为观测角度向量。
优选地，估计矢量
Figure BDA0002202810440000054
分别根据公式(7)确定：
其中，
Figure BDA0002202810440000056
表示观测平台与目标星间的距离，因距离大小不影响矢量方向，计算中，
Figure BDA0002202810440000057
可取任意正整数；
Figure BDA0002202810440000058
对应估计的目标星俯仰角和方位角。
优选地，根据公式(8)确定，其推算过程如下：
当t0时刻估计初值为
Figure BDA0002202810440000061
时，可以得到ti时刻目标星在观测平台坐标系下的估计位置
Figure BDA0002202810440000063
继而由式(8)得到估计角度向量
Figure BDA0002202810440000064
Figure BDA0002202810440000065
其中，
Figure BDA0002202810440000066
和为目标在观测平台坐标系下的估计空间位置。
优选地，所述步骤S7中的vk,g根据变异公式(9)确定：
vk,g＝Xbest,g+F·(Xr1,g-Xr2,g) (9)
其中，F为变异缩放因子，其范围为(0,1)；r1,r2是0到NP-1之间的随机正整数,r1≠r2≠k。
优选地，所述步骤S8中的uk,g根据交叉公式(10)确定：
Figure BDA0002202810440000068
其中jrand为0到NP-1间的随机整数，Cr为变异率。
优选地，所述步骤S9中的xk,g+1根据公式(11)确定：
Figure BDA0002202810440000069
与现有技术相比，本发明的有益效果是：
(1)本发明提供的评价函数可根据任务需求灵活设计，不增加求解难度。
(2)本发明最小二乘批处理精定轨的定轨思路，设计了基于差分进化的仅测角精定轨方法。差分进化定轨相较现有的轨道改进方法在目标状态解算环节进行了创新，使用差分进化方法直接在解空间搜索时估计残差方差全局最小的目标状态，省略了初定轨需求，避免了运算中的线性化误差，并且不需要对正规矩阵进行计算，避免了计算过程引入的误差，可以提高定轨精度。
(3)本发明针对天基观测系统存在的问题，基于轨道动力学方程，提出了一种基于差分进化的定轨方法，该方法利用卫星仅测角信息，可以对空间非合作进行精密定轨。相较于传统最小二乘批处理精密定轨方法，该算法不需要初定轨信息，并避免了求解过程中的发散问题及线性化引入误差。仿真验证结果表明，当存在观测噪声时，增加观测弧长可使定轨结果更加接近真实值，从而达到精密定轨的目标。
附图说明
图1为本发明提供的基于差分进化的仅测角定轨流程图；
图2本发明提供的观测平台轨道坐标系下测角示意图；
图3为本发明提供的地心赤道坐标系及观测平台轨道坐标系示意图；
图4为地心坐标系下的卫星示意图；
图5为观测平台轨道坐标系下测角示意图。
具体实施方式
下面结合附图针对本发明作进一步实例描述：
在实际的轨道计算中，常用的方法就是根据已知的用初始向量和卫星动力学方程进行计算来得到卫星定轨信息，本发明通过直接计算初始向量，符合实际轨道计算的需求。此外，本发明提供的不需要计算传统方法中的正规矩阵及矩阵求逆，不受中观测向量和预测向量长度的限制，这样可以充分利用观测向量信息来提高定轨精度，而且评价函数中直接使用实际的轨道预测信息，使得预测结果更能满足实际需求。同时为克服寻优空间维数升高后，传统的遗传算法计算效率不高的问题，本发明采取差分进化算法来进行寻优，且根据定轨条件对寻优空间范围进行预处理，进一步加快收敛速度。
本发明具体采用下述技术方案：采用下述技术方案：一种基于差分进化的仅测角定轨方法，其特征在于，包括如下步骤：
S1：将卫星差t0时刻的位置速度状态作为参数，记为
Figure BDA0002202810440000081
并设置差分化进化算法参数，所述参数包括变量维数D、种群个数NP、变异缩放因子F；变异率Cr；最大迭代代数gmax，其中，D为6；
S2：根据步骤S1中设置的变量维数D、种群个数NP进行种群初始化，生成初始种群，种群个体的目标星轨道参数预测值为Xk,g，k＝1…NP；
S3：对步骤S2生成的初始种群，计算以种群个体Xk,g为初始状态时，在各观测时刻观测方向与估计方向夹角θi(i＝1,…,n，n为观测点个数)，将夹角θi作为i时刻的估计误差，将各时刻误差序列表示为向量形式：r＝[θ0,θ1,…,θn]T。
S4：通过最小二乘法得到估计值残差定义，设置估计值残差为RSM，通过式(1)对步骤S2生成种群的每个个体进行残差值计算；
RMS＝rTr (1)
其中r＝[θ0,θ1,…,θn]T。
S5：检查并剔除步骤S4中误差向量r中的异常值，若误差向量r的第i项|ri|≥3σ，其中，即看作异常值；
其中，σ2为r各分量平方的均值，
S6：使用步骤S4得到的残差值RMS和步骤S5得到的异常值计算评价函数，得到Xk,g的评价函数值ffitness(X)：评价函数如式(2)所示，
其中，m为异常值个数，n-m为非异常值个数，routlier为差值向量r中m个异常值组成的向量；
S7：变异，通过比较步骤S6中获得的Xk,g的评价函数值ffitness(X)，得出评价函数最优状态向量Xbest,g，并通过最优状态向量Xbest,g进行变异得到交叉变异状态向量vk,g；
S8：交叉，通过交叉变异状态向量vk,g和原状态向量Xk,g中不同位置元素，产生新的候选状态向量uk,g；
S9：选择，通过对比候选状态向量uk,g和原状态向量Xk,g的评价函数，生成下一代种群的状态向量Xk,g+1；
S10：终止条件，设第g代最佳状态向量Xbest,g对应的评价函数为fbest,g，设最大进化代数为gmax；当达到以下条件之一时，差分进化计算过程终止：
(1)fbest,g-fbest,g-1＜δ，δ正整数
(2)g＞gmax
取终止计算时的Xbest,g作为定轨结果。
其中，本发明提供的步骤S2中种群的初始化过程具体为：
设有效观测距离为d，观测平台轨道参数为
Figure BDA0002202810440000091
在状态空间中，rmin表示航天器与地心距离的最小值，根据rmin和目标偏心率范围可以计算出速度最大值vmax，即坐标轴各方向的速度取值范围为[-vmax,vmax]；
根据公式(3)进行种群初始化过程，
Figure BDA0002202810440000101
其中，状态向量上界bU＝(xT+d,yT+d,zT+d,vmax,vmax,vmax)T，状态变量下界bL＝(xT-d,yT-d,zT-d,-vmax,-vmax,-vmax)T，Xj,k,g为第g代第k个个体的第j项；P为种群，表示所有的状态向量，g为当前迭代代数；D为变量维数。
其中，参图2所示，图2为观测平台轨道坐标系下测角示意图，步骤S3中的观测方向与估计方向夹角θi通过公式(4)确定：
Figure BDA0002202810440000102
其中，设ti时刻观测平台坐标系下观测矢量为Li，估计矢量为
Figure BDA0002202810440000103
其中，观测矢量L根据公式(5)确定：
L＝ρ[cosαcosβ,sinαcosβ,sinβ]T,ρ＞0 (5)
其中，ρ表示观测平台与目标星间的距离，因距离大小不影响矢量方向，计算中，ρ可取任意正整数；α、β为实际观测的俯仰角和方位角。
α、β根据公式(6)确定，其推算过程如下：
当t0时刻估计初值为X0时，可以得到ti时刻目标星在观测平台坐标系下的位置rs＝[xs,ys,zs]T时，对应的观测角的可通过式(6)计算：
Figure BDA0002202810440000104
其中，xs,ys和zs为目标在观测平台坐标系下的空间位置，其中，H＝[αβ]T为观测角度向量。
其中，估计矢量
Figure BDA0002202810440000105
分别根据公式(7)确定：
Figure BDA0002202810440000111
其中，
Figure BDA0002202810440000112
表示观测平台与目标星间的距离，因距离大小不影响矢量方向，计算中，
Figure BDA0002202810440000113
可取任意正整数；
Figure BDA0002202810440000114
对应估计的目标星俯仰角和方位角。
Figure BDA0002202810440000115
根据公式(8)确定，其推算过程如下：
当t0时刻估计初值为
Figure BDA0002202810440000116
时，可以得到ti时刻目标星在观测平台坐标系下的估计位置
Figure BDA0002202810440000117
Figure BDA0002202810440000118
继而由式(8)得到估计角度向量
Figure BDA00022028104400001110
其中，
Figure BDA00022028104400001111
和
Figure BDA00022028104400001112
为目标在观测平台坐标系下的估计空间位置。
其中，本发明提供的所述步骤S7中的vk,g根据变异公式(9)确定：
vk,g＝Xbest,g+F·(Xr1,g-Xr2,g) (9)
其中，F为变异缩放因子，其范围为(0,1)；r1,r2是0到NP-1之间的随机正整数,r1≠r2≠k。
本发明提供的所述步骤S8中的uk,g根据交叉公式(10)确定：
Figure BDA00022028104400001113
其中jrand为0到NP-1间的随机整数，Cr为变异率。
本发明提供的所述步骤S9中的xk,g+1根据公式(11)确定：
其中，步骤S3对步骤S2生成的初始种群，计算以种群个体Xk,g为初始状态时，在各观测时刻观测方向与估计方向夹角θi(i＝1,…,n，n为观测点个数)，将夹角θi作为i时刻的估计误差，其中，本发明将Xk,g作为t0时刻初值，得到估计轨道，继而由估计轨道得到估计的目标在观测平台坐标系中位置，得到估计位置矢量。计算估计位置矢量
Figure BDA0002202810440000121
和实际观测位置矢量Li两个矢量各时刻夹角θi，得到定义的误差序列；
其中，本发明估计位置矢量和实际观测位置矢量Li需要先计算α、β及
Figure BDA0002202810440000123
Figure BDA0002202810440000124
α、β及
Figure BDA0002202810440000125
分别根据公式(6)和(8)确定，具体的推算过程如下：
参图3所示，图3为地心赤道坐标系及观测平台轨道坐标系示意图。地心赤道惯性坐标系OEXEYEZE的原点在地心OE，XE轴在赤道面内指向春分点Υ，YE轴在赤道平面内与XE轴垂直指向动，ZE轴与XE轴、YE轴构成右手正交坐标系。观测平台轨道坐标系OCXCYCZC的原点OC为观测平台在地心赤道坐标系中的空间位置，XC轴沿地心OE指向观测平台OC，YC轴在轨道平面上与XC轴垂直，指向航天器的速度方向，ZC轴与XC轴、YC轴组成右手正交坐标系，垂直于观测平台的轨道平面。
将卫星在地心坐标系中的位置矢量用r表示，r＝[x,y,z]T。卫星受地球引力及其他天体带来的摄动力，根据牛顿定律，可以得到卫星的运动方程：
Figure BDA0002202810440000126
式(12)中，
Figure BDA0002202810440000127
为卫星到地心距离，G为万有引力常数，M为地球质量，m为卫星质量，F为摄动力。将式(12)表示为加速度形式：
Figure BDA0002202810440000128
式(13)中，μ＝GM为地心引力常数，
Figure BDA0002202810440000131
为地球引力外其他天体摄动力带来的加速度。记在地心坐标系中，方程(13)可沿坐标轴分解为：
Figure BDA0002202810440000133
其中，ax，ay，az分别为a的在X轴、Y轴、Z轴的分量。
将卫星在地心坐标系的速度位置信息表示为使用式(4)构造X的微分方程，如式(14)所示：
Figure BDA0002202810440000135
根据式(15)，可以得到X的状态方程，如式(16)所示：
Figure BDA0002202810440000136
给定t0时刻卫星位置速度X0时，由状态方程(16)，通过数值积分运算，可以确定未来任意时刻的卫星的位置速度X(t)。
观测平台SC与目标星St在地心坐标系中的关系示意图如图4所示，图4为地心坐标系下的卫星示意图，图中rc，rt分别为观测平台和目标星的位置矢量，LI为观测平台到目标星的相对位置矢量。
如前所述，当t0时刻目标星在地心坐标系下空间位置速度为Xt(t0)时，可以求得t0后任意时刻目标星在地心坐标系下的空间位置rt。在地心坐标系中，目标星与观测平台相对位置为
LI＝rt-rc (17)
在仅测角定轨任务中，观测的方位角α和俯仰角β可通过目标星在观测平台轨道坐标系下的位置矢量计算。参图4所示，图4为位置矢量L与测角关系示意图。L可通过地心坐标系下目标星与观测平台的相对位置矢量LI经过坐标转换计算：
L＝R·LI (18)
式(18)中，R为地心坐标系到观测平台坐标系的坐标转换矩阵。坐标转换矩阵的求解过程如下：
地心坐标系下，设i，j，k分别为XE轴，YE轴和ZE轴方向的单位向量，其中，i＝[1,0,0]T，j＝[0,1,0]T，k＝[0,0,1]T，地心坐标系下任意位置矢量r＝[x,y,z]T可表示为
Figure BDA0002202810440000141
设ic，jc，kc分别为观测卫星轨道坐标系XC轴，YC轴，ZC轴单位向量在地心坐标系下的向量表示，由图2观测星轨道坐标系各坐标轴的定义，可以得到：
Figure BDA0002202810440000142
式(20)中，rt为地心坐标系下观测星的位置矢量。设地心坐标系下位置矢量r＝[x,y,z]T在观测卫星轨道坐标系下的坐标为[xc,yc,zc]T，可以得到如下关系：
Figure BDA0002202810440000151
由式(21)可以直接得到任意位置矢量从地心坐标系到观测平台轨道坐标系的坐标变换，如式(22)所示：
Figure BDA0002202810440000152
由于[ic,jc,kc]为单位正交矩阵，[i,j,k]为单位矩阵，[ic,jc,kc]-1[i,j,k]＝[ic,jc,kc]T，将[ic,jc,kc]T简写为R，式(22)可表示为
Figure BDA0002202810440000153
由式(23)可知，R即为所求地心坐标系到观测平台坐标系的坐标转换矩阵，其求解如式(20)所示。
参图2所示，图2为观测平台轨道坐标系下测角示意图
当某时刻目标星在观测平台坐标系下位置为rs＝[xs,ys,zs]T时，对应的观测角的可通过(6)计算：
Figure BDA0002202810440000154
其中，H＝[αβ]T为观测角度向量。
根据式(16)、(17)和(20)，当t0时刻估计初值为
Figure BDA0002202810440000161
时，可以得到ti时刻目标星在观测平台坐标系下的估计位置
Figure BDA0002202810440000162
继而由式(8)得到估计角度向量
Figure BDA0002202810440000163
如图2所示。
实施例1
仿真验证-定轨实验
本发明选取FengYun1C空间碎片作为观测目标，其轨道数据取自公开的TLE数据集[celestrak.com]。FengYun1C的轨道高度与太阳同步卫星轨道接近，观测平台轨道参照风云一号选取，观测平台半长轴a＝7237km，偏心率为e＝0，轨道倾角98.88度。实验随机选取了空间分布不同的空间碎片作为观测目标。
实验中轨道数据考虑了J2摄动影响，使用STK生成卫星轨道位置、速度数据和角度数据，并在matlab中进行了方法验证。角度数据观测精度为5”，采样频率1HZ，有效观测距离5000km，实验选取连续观测时间不小于10min的角度序列作为观测序列。试验中噪声分布为标准正态分布。
本实施例中，以进化过程最快收敛目标，经过测试对比，最终选取参数如下：
种群个数NP＝50；变异缩放因子F＝0.9；变异率Cr＝1；最大迭代代数g＝400。(1)定轨精度与相对距离关系分析
实验将不同空间分布的目标卫星定轨结果进行了比较，表1为观测弧长10min时，定轨结果与真实值的差值情况。
表1卫星定轨误差统计表
Figure BDA0002202810440000165
Figure BDA0002202810440000171
因观测误差的分布具有随机性，单次定轨结果往往不能反映实际定轨效果，为保证结果的可靠性，试验对各卫星取50次定轨误差的平均值作为结果，每次计算重新生成观测误差。
(2)定轨精度与观测弧长关系
观测序列较短时，观测噪声一般不符合正态分布，使用噪声干扰的观测数据进行估计，结果通常会偏离真实值。Vallado的研究表明，通过增加观测弧长，可以提高定轨精度，并提出实现精密定轨使用的观测弧长应不小于5min[1998]，但弧长的增加会增加计算时间，因此实验对观测弧长和定轨精度的关系进行了研究。表2、3显示了各卫星在观测弧长分别为300s、400、500s时的定轨结果统计，显然较长的观测弧段能提供更好的定轨精度。
表2目标星1定轨精度与观测弧长关系
Figure BDA0002202810440000172
表3目标星2定轨精度与观测弧长关系
(3)轨道预测
使用观测时长10min的定轨结果，通过动力学方程计算预测目标1h后的状态，位置及速度差值如表4所示。
表4轨道估计误差统计表
Figure BDA0002202810440000182
仿真验证结果表明，利用本发明提供的定轨方法进行定轨，可以达到精密定轨的目标。且当存在观测噪声时，增加观测弧长可使定轨结果更加接近真实值，从而达到精密定轨的目标。
以上所述仅为本发明的实施例子，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (10)

1.一种天基仅测角差分进化一次精确定轨方法，其特征在于，包括如下步骤：
S1：将卫星差t0时刻的位置速度状态作为参数，记为并设置差分化进化算法参数，所述参数包括变量维数D、种群个数NP、变异缩放因子F；变异率Cr；最大迭代代数gmax，其中，D为6；
S2：根据步骤S1中设置的变量维数D、种群个数NP进行种群初始化，生成初始种群，种群个体的目标星轨道参数预测值为Xk,g，k＝1…NP；
S3：对步骤S2生成的初始种群，计算以种群个体Xk,g为初始状态时，在各观测时刻观测方向与估计方向夹角θi(i＝1,…,n，n为观测点个数)，将夹角θi作为i时刻的估计误差，将各时刻误差序列表示为向量形式：r＝[θ0,θ1,…,θn]T；
S4：通过最小二乘法得到估计值残差定义，设置估计值残差为RSM，通过式(1)对步骤S2生成种群的每个个体进行残差值计算；
RMS＝rTr (1)
其中r＝[θ0,θ1,…,θn]T；
S5：检查并剔除步骤S4中误差向量r中的异常值，若误差向量r的第i项|ri|≥3σ，其中，即看作异常值；
其中，σ2为r各分量平方的均值，
Figure FDA0002202810430000012
S6：使用步骤S4得到的残差值RMS和步骤S5得到的异常值计算评价函数，得到Xk,g的评价函数值ffitness(X)：评价函数如式(2)所示，
Figure FDA0002202810430000013
其中，m为异常值个数，n-m为非异常值个数，routlier为差值向量r中m个异常值组成的向量；
S7：变异，通过比较步骤S6中获得的Xk,g的评价函数值ffitness(X)，得出评价函数最优状态向量Xbest,g，并通过最优状态向量Xbest,g进行变异得到交叉变异状态向量vk,g；
S8：交叉，通过交叉变异状态向量vk,g和原状态向量Xk,g中不同位置元素，产生新的候选状态向量uk,g；
S9：选择，通过对比候选状态向量uk,g和原状态向量Xk,g的评价函数，生成下一代种群的状态向量Xk,g+1；
S10：终止条件，设第g代最佳状态向量Xbest,g对应的评价函数为fbest,g，设最大进化代数为gmax；当达到以下条件之一时，差分进化计算过程终止：
(1)fbest,g-fbest,g-1＜δ，δ正整数
(2)g＞gmax
取终止计算时的Xbest,g作为定轨结果。
2.如权利要求1所述的定轨方法，其特征在于，所述步骤S2中种群的初始化过程具体为：
设有效观测距离为d，观测平台轨道参数为
Figure FDA0002202810430000021
在状态空间中，rmin表示航天器与地心距离的最小值，根据rmin和目标偏心率范围可以计算出速度最大值vmax，即坐标轴各方向的速度取值范围为[-vmax,vmax]；
根据公式(3)进行种群初始化过程，
其中，状态向量上界bU＝(xT+d,yT+d,zT+d,vmax,vmax,vmax)T，状态变量下界bL＝(xT-d,yT-d,zT-d,-vmax,-vmax,-vmax)T，Xj,k,g为第g代第k个个体的第j项；P为种群，表示所有的状态向量，g为当前迭代代数；D为变量维数。
3.如权利要求1所述的定轨方法，其特征在于，步骤S3中的观测方向与估计方向夹角θi通过公式(4)确定：
其中，设ti时刻观测平台坐标系下观测矢量为Li，估计矢量为
Figure FDA0002202810430000033
4.如权利要求3所述的定轨方法，其特征在于，观测矢量Li根据公式(5)：
L＝ρ[cosαcosβ,sinαcosβ,sinβ]T,ρ＞0 (5)
其中，ρ表示观测平台与目标星间的距离，因距离大小不影响矢量方向，计算中，ρ可取任意正整数；α、β为实际观测的俯仰角和方位角。
5.如权利要求4所述的定轨方法，其特征在于，α、β根据公式(6)确定，其推算过程如下：
当t0时刻估计初值为X0时，可以得到ti时刻目标星在观测平台坐标系下的位置rs＝[xs,ys,zs]T时，对应的观测角的可通过式(6)计算：
Figure FDA0002202810430000034
其中，xs,ys和zs为目标在观测平台坐标系下的空间位置，其中，H＝[α β]T为观测角度向量。
6.如权利要求3所述的定轨方法，其特征在于，估计矢量
Figure FDA0002202810430000041
分别根据公式(7)确定：
Figure FDA0002202810430000042
其中，
Figure FDA0002202810430000043
表示观测平台与目标星间的距离，因距离大小不影响矢量方向，计算中，
Figure FDA0002202810430000044
可取任意正整数；
Figure FDA0002202810430000045
对应估计的目标星俯仰角和方位角。
7.如权利要求4所述的定轨方法，其特征在于，
Figure FDA0002202810430000046
根据公式(8)确定，其推算过程如下：
当t0时刻估计初值为
Figure FDA0002202810430000047
时，可以得到ti时刻目标星在观测平台坐标系下的估计位置
Figure FDA0002202810430000048
Figure FDA0002202810430000049
继而由式(8)得到估计角度向量
Figure FDA00022028104300000410
Figure FDA00022028104300000411
其中，
Figure FDA00022028104300000412
和为目标在观测平台坐标系下的估计空间位置。
8.如权利要求1所述的定轨方法，其特征在于，所述步骤S7中的vk,g根据变异公式(9)确定：
vk,g＝Xbest,g+F·(Xr1,g-Xr2,g) (9)
其中，F为变异缩放因子，其范围为(0,1)；r1,r2是0到NP-1之间的随机正整数,r1≠r2≠k。
9.如权利要求1所述的定轨方法，其特征在于，所述步骤S8中的uk,g根据交叉公式(10)确定：
Figure FDA0002202810430000051
其中jrand为0到NP-1间的随机整数，Cr为变异率。
10.如权利要求1所述的定轨方法，其特征在于，所述步骤S9中的xk,g+1根据公式(11)确定：
Figure FDA0002202810430000052
CN201910870972.9A 2019-09-16 2019-09-16 一种天基仅测角差分进化一次精确定轨方法 Active CN110779531B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910870972.9A CN110779531B (zh) 2019-09-16 2019-09-16 一种天基仅测角差分进化一次精确定轨方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910870972.9A CN110779531B (zh) 2019-09-16 2019-09-16 一种天基仅测角差分进化一次精确定轨方法

Publications (2)

Publication Number Publication Date
CN110779531A true CN110779531A (zh) 2020-02-11
CN110779531B CN110779531B (zh) 2020-11-20

Family

ID=69383486

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910870972.9A Active CN110779531B (zh) 2019-09-16 2019-09-16 一种天基仅测角差分进化一次精确定轨方法

Country Status (1)

Country Link
CN (1) CN110779531B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111637895A (zh) * 2020-05-11 2020-09-08 北京控制工程研究所 一种基于q学习的导航观测目标选取方法
CN115790607A (zh) * 2023-01-31 2023-03-14 南京航空航天大学 一种基于短弧历史数据的非合作目标机动特征检测方法
CN116258090A (zh) * 2023-05-16 2023-06-13 中国地质大学（武汉） 基于双阶段信息迁移的差分演化深空轨道设计方法及系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160047878A1 (en) * 2011-10-25 2016-02-18 The Boeing Company Combined orbit and attitude determination system and methods
CN108287334A (zh) * 2018-02-06 2018-07-17 西安四方星途测控技术有限公司 一种基于rcs测量数据的自旋卫星姿态估计方法及系统
CN109284579A (zh) * 2018-11-12 2019-01-29 东南大学 基于改进差分进化算法的非线性系统参数辨识方法
CN109631911A (zh) * 2018-12-17 2019-04-16 浙江大学 一种基于深度学习目标识别算法的卫星姿态转动信息确定方法
CN109885872A (zh) * 2019-01-10 2019-06-14 杭州电子科技大学 一种基于差分进化算法的均匀面阵稀疏优化方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160047878A1 (en) * 2011-10-25 2016-02-18 The Boeing Company Combined orbit and attitude determination system and methods
CN108287334A (zh) * 2018-02-06 2018-07-17 西安四方星途测控技术有限公司 一种基于rcs测量数据的自旋卫星姿态估计方法及系统
CN109284579A (zh) * 2018-11-12 2019-01-29 东南大学 基于改进差分进化算法的非线性系统参数辨识方法
CN109631911A (zh) * 2018-12-17 2019-04-16 浙江大学 一种基于深度学习目标识别算法的卫星姿态转动信息确定方法
CN109885872A (zh) * 2019-01-10 2019-06-14 杭州电子科技大学 一种基于差分进化算法的均匀面阵稀疏优化方法

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111637895A (zh) * 2020-05-11 2020-09-08 北京控制工程研究所 一种基于q学习的导航观测目标选取方法
CN111637895B (zh) * 2020-05-11 2021-10-01 北京控制工程研究所 一种基于q学习的导航观测目标选取方法
CN115790607A (zh) * 2023-01-31 2023-03-14 南京航空航天大学 一种基于短弧历史数据的非合作目标机动特征检测方法
CN115790607B (zh) * 2023-01-31 2023-05-12 南京航空航天大学 一种基于短弧历史数据的非合作目标机动特征检测方法
CN116258090A (zh) * 2023-05-16 2023-06-13 中国地质大学（武汉） 基于双阶段信息迁移的差分演化深空轨道设计方法及系统
CN116258090B (zh) * 2023-05-16 2023-08-18 中国地质大学（武汉） 基于双阶段信息迁移的差分演化深空轨道设计方法及系统

Also Published As

Publication number Publication date
CN110779531B (zh) 2020-11-20

Similar Documents

Publication Publication Date Title
CN110779531B (zh) 一种天基仅测角差分进化一次精确定轨方法
US8639476B2 (en) Process for estimation of ballistic missile boost state
CN107300700B (zh) 敏捷合成孔径雷达卫星聚束模式姿态机动需求计算方法
CN108562295A (zh) 一种基于同步卫星二体模型的三站时差定轨方法
Gong et al. Maneuver-free approach to range-only initial relative orbit determination for spacecraft proximity operations
CN105701283B (zh) 地球非球形摄动作用下自由段弹道误差传播的分析方法
Liu et al. Analytical investigations of quasi-circular frozen orbits in the Martian gravity field
CN105825058B (zh) 超稀疏雷达数据摄动补偿初轨计算方法
CN113587926A (zh) 一种航天器空间自主交会对接相对导航方法
CN115143955B (zh) 基于天文测角数据确定地球同步轨道带航天器初轨的方法
Verma et al. Prospects of dynamical determination of General Relativity parameter beta and solar quadrupole moment J2 with asteroid radar astronomy
Chen et al. Optimal orbit design for repeat-pass multi-baseline observation of tomographic SAR satellites
Zhu et al. Mean-shift clustering approach to the tracklets association with angular measurements of resident space objects
CN115774928B (zh) 基于改进Laplace模型的空间碎片短弧仅测角初定轨优化方法
Cordelli et al. Fusion of laser ranges and angular measurements data for LEO and GEO space debris orbit determination
Liu et al. Applying Lambert problem to association of radar-measured orbit tracks of space objects
Reihs et al. Tracklet-based correlation of combined radar and optical measurements
Davis et al. Identifying near-term missions and impact keyholes for asteroid 99942 apophis
Tartakovsky et al. Efficient Estimation and Decision-Making Methods for Short Track Identification and Orbit Determination
Zhao et al. A Fast Algorithm for Determining the Optimal Navigation Star for Responsive Launch Vehicles
Pardal et al. Comparing the extended and the sigma point Kalman filters for orbit determination modeling using GPS measurements
Wu et al. Random Error Estimation and Propagation Analysis for Satellites’ Initial Positions
CN113282096B (zh) 地球静止轨道博弈航天器相对位置非线性误差的控制方法
Hill et al. Liaison navigation in the sun-earth-moon four-body problem
Bae et al. The GLAS Algorithm Theoretical Basis Document for Precision Attitude Determination (PAD)

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