技术领域
[0001] 本
发明涉及
航天器编队飞行空间测量领域，特别是涉及基于星敏感器和星间链路的HEO卫星编队飞行自主导航方法。
背景技术
[0002] 大椭圆轨道(High Elliptical Orbit，HEO)卫星的轨道偏心率大于0.6，近地点高度在300千米～1000千米之间，远地点处的高度要高于地球静止轨道卫星(36000千米)，其原理图如图1-1所示。HEO卫星具有
覆盖区域广，覆盖时间长的特点，因此在很多领域都得到了应用，例如导弹预警，气象探测，导航
定位，空间科学探索等。早期的HEO卫星主要应用于军事领域，其中的典型代表有俄罗斯的“闪电”系列卫星；美国的“折叠椅”，“军号”卫星。可以说，尽管目前世界各国研究的重点仍然是低轨卫星和中轨卫星，但HEO卫星由于其自身独特的轨道特点，可以做为现有卫星的补充，必将拥有广泛的应用前景。
[0003] HEO卫星群在空间运行过程中所经历的空间环境比较复杂，要想实现其高
精度的自主导航，必须解决两方面问题：一方面，当噪声的分布并不确定时，所采用的滤波
算法是否具有较强的自适应能
力使得滤波结果收敛，并且保证足够高的精度；另一方面，HEO卫星的速度和
加速度在近地点附近会出现较为剧烈的变化，即存在着高速度高动态的问题，它会导致滤波结果出现
波动的情况。
[0004] 为了HEO卫星进行编队飞行以及共位控制，必须首先能实时获得卫星的
位置和
姿态信息，并且不能对临星产生干扰，由于HEO卫星远地点通常在36000km高度，HEO卫星存在导航观测
信号不足的问题：(1)常用GNSS导航方式存在导航信号弱、地球遮挡和可见卫星少等严重问题；(2)其它导航方式：地
磁场无法使用，雷达高度计仅限于近地轨道导航，都难以作为观测信息满足整个轨道周期导航要求，因此这就使得研究HEO卫星引入新观测方法显得迫切重要。
发明内容
[0005] 发明目的：本发明针对HEO卫星编队飞行观测信息不足导致导航精度较低的问题，提出一种基于星敏感器和星间链路的HEO卫星编队飞行自主导航方法，为编队飞行的HEO卫星提供高精度相对观测信息。
[0006] 技术方案：为达到此目的，本发明采用以下技术方案：
[0007] 本发明所述的基于星敏感器和星间链路的HEO卫星编队飞行自主导航方法，包括以下步骤：
[0008] S1：以HEO卫星对地观测为任务需求，将两个HEO卫星分别设为主星和子星，设计主星和子星编队飞行构型及轨道参数；
[0009] S2：根据地心惯性
坐标系下卫星相对轨道动力学模型，建立自主
导航系统状态模型；
[0010] S3：根据步骤S1计算主星和子星相对距离，判断子星是否满足星敏感器观测距离要求：如果满足，则进入步骤S4；否则，进入步骤S12；
[0011] S4：根据步骤S1解算太阳、地球和子星三者的位置关系，判断子星是否处在太阳光照区：如果是，则进入步骤S5；否则，进入步骤S12；
[0012] S5：根据步骤S1解算地球、主星和子星三者的位置关系，判断地球否进入星敏感器视场：如果是，则进入步骤S6；否则，进入步骤S12；
[0013] S6：根据步骤S1计算子星可视星等，判断子星可视星等是否小于星敏感器可观测
阈值：如果是，则进入步骤S7；否则，进入步骤S12；
[0014] S7：根据步骤S1计算主星相对子星方向矢量与星敏感器光轴指向夹
角，判断子星是否在星敏感器视场范围内：如果是，则进入步骤S8；否则，通过万向轴调整星敏感器光轴指向，再次按照步骤S7判断子星是否在星敏感器视场范围内，是则进入步骤S8，否则进入S12；
[0015] S8：根据步骤S7计算子星在星敏感器二维像面阵坐标，判断子星是否在星敏感器二维像面阵内：如果是，则进入步骤S9；否则，进入步骤S12；
[0016] S9：计算主星相对子星的理论方向矢量、方位角以及
俯仰角；
[0017] S10：调整星敏感器真实光轴指向与理论方向矢量一致，对子星进行真实观测，计算子星相对卫星真实方向矢量，建立以单位方向矢量和距离为观测量的观测方程；
[0018] S11：对步骤S2所建立的状态模型和步骤S10所建立的观测方程离散化，利用Unscented卡尔曼滤波算法估计卫星位置和速度；
[0019] S12：结束观测。
[0020] 进一步，所述步骤S2中，建立自主导航系统状态模型的过程如下：
[0021] 在地心惯性坐标系下，当主星位置距离大于子星与主星相对距离的时候，建立卫星相对目标子星轨道动力学模型：
[0022]
[0023]
[0024] 其中，δr(10)为主星相对子星方向矢量，δv(10)为主星相对子星速度矢量， 为卫星位置随时间的一阶导数， 为卫星速度随时间的一阶导数，r(0)为主星位置矢量，μe为地球引力常数，af为摄动加速度；
[0025] 定义状态变量x＝[(δr(10))T(δv(10))T]T，建立自主导航系统状态模型；
[0026]
[0027] 其中，f为系统非线性连续状态转移函数，wt为状态噪声。
[0028] 进一步，所述步骤S3中，判断子星否满足星敏感器观测距离要求的过程如下：计算主星相对子星距离δr(10)，判断δr(10)是否满足式(3)所示的条件：
[0029] Lmin≤δr(10)≤Lmax (3)
[0030] 其中，δr(10)＝|δr(10)|＝|r(1)-r(0)|，r(0)为主星位置矢量，r(1)为子星位置矢量，Lmin星间观测所需最小距离，Lmax为星间观测所需最大距离。
[0031] 进一步，所述步骤S4中，判断子星是否处在太阳光照区的过程如下：分析地球阴影范围以及子星运行穿过该阴影范围的临界条件，设子星位置矢量r(1)与太阳位置矢量r(sun)夹角为ψ，子星进入和离开地球阴影范围的临界夹角为ψcri，则子星处在太阳光照区需要满足条件：
[0032] ψ＜ψcri (4)。
[0033] 进一步，所述步骤S5中，判断地球是否进入星敏感器视场的过程如下：设主星位置矢量r(0)和主星相对子星方向矢量δr(10)的夹角为θ，被地球遮挡导致背景光线过弱的临界条件是主星相对子星方向矢量δr(10)与地球边缘相切，定义此临界夹角为θcri，则地球未进入星敏感器视场条件为：
[0034] θ＞θcri (5)。
[0035] 进一步，所述步骤S6中，子星可视星等根据式(6)计算得出：
[0036]
[0037] 式(6)中，m为子星可视星等，M为子星绝对星等，根据式(7)计算得到；|r(sun0)|是太阳与子星之间的距离；ξ是主星相对子星方向矢量δr(10)与太阳相对子星方向矢量r(sun1)的夹角，p(ξ)是
相位积分，根据式(8)计算得到；d0是地球与太阳之间的平均距离；
[0038]
[0039] 式(7)中，msun是太阳的可视星等，rd为被观测天体的半径，a是被观测天体的反射率；
[0040]
[0041] 进一步，所述步骤S7中，判断子星是否在星敏感器视场范围内的过程如下：如果则子星在星敏感器视场范围内；否则，则子星不在星敏感器视场范围内；如式(9)所示，FOV为星敏感器视场角；
[0042]
[0043] 式(9)中，δr(10)为主星相对子星方向矢量， 为星敏感器方向矢量， 是地心惯性坐标系相对本体坐标系姿态转换矩阵；
[0044] 所述步骤S8中，判断子星是否在星敏感器二维像面阵内的过程如下：如果满足式(10a)和(10b)，则子星在星敏感器二维像面阵内；
[0045]
[0046]
[0047] 其中， 为主星相对子星方向矢量δr(10)投影在星敏感器二维像面阵的坐标，IPlongth为二维像面阵的长度，IPwidth为二维像面阵的宽度。
[0048] 进一步，所述步骤S9中，主星相对子星的理论方向矢量、方位角和俯仰角分别根据式(11)、(12a)和(12b)计算得到：
[0049]
[0050] 式(11)中，δr(10)为主星相对子星的理论方向矢量，也称为主星相对子星方向矢量； 为主星相对子星单位方向矢量；
[0051]
[0052]
[0053] 其中， 是地心惯性坐标系相对本体坐标系姿态转换矩阵，α为方位角，δ为俯仰角。
[0054] 进一步，所述步骤S10中，建立的观测方程为：
[0055]
[0056] 式(13)中， 为主星相对子星单位方向矢量的真实测量值， 为主星相对子星距离的真实测量值。
[0057] 进一步，所述步骤S11的具体过程如下：
[0058] 根据式(14a)对步骤S2所建立的状态模型进行离散化：
[0059]
[0060] 式(14a)中，xk∈RL为第k个状态变量，xk+1为第k+1个状态变量，wk∈N(0,Qk)为第k个过程噪声，Qk为系统噪声强度，f为系统非线性连续状态转移函数；
[0061] 根据式(14b)对步骤S10所建立的观测方程进行离散化：
[0062] yk＝g(xk)+vk (14b)
[0063] 式(14b)中，yk∈RM为第k个输出向量，g为观测方程，vk∈N(0,Rk)为第k个测量噪声，Rk为观测噪声强度，wk和vk不相关；
[0064] 利用Unscented卡尔曼滤波算法估计卫星位置和速度，步骤如下：
[0065] S11.1：根据式(15a)、(15b)和(15c)进行Unscented变换：
[0066]
[0067]
[0068]
[0069] 其中，为xk的均值， 为xk的方差，λ＝α2(1+κ)-1为标量参数，α决定了sigma点在周围的分布，κ为标量参数，L为sigma点维数，χi,k-1是中间变量；
[0070] S11.2：预测过程：根据式(16a)、(16b)、(16c)、(16d)、(16e)和(16f)进行预测；
[0071] χi,k/k-1＝f(χi,k-1) (16a)
[0072]
[0073]
[0074]
[0075]
[0076]
[0077] 其中，Wi和 分别为U变换计算均值和方差所用的权系数，χi,k/k-1、Yi,k/k-1和 都是中间变量；
[0078] S11.3：更新过程：根据式(17a)、(17b)、(17c)、(17d)和(17e)进行更新；
[0079]
[0080]
[0081]
[0082]
[0083]
[0084] 其中， 和Kk都是中间变量；
[0085] S11.4：返回步骤S11.1进行下一个周期的滤波。
[0086] 有益效果：本发明公开了一种基于星敏感器和星间链路的HEO卫星编队飞行自主导航方法，有效解决卫星编队飞行观测信息不足所导致的导航精度较低的问题，相对于
现有技术，具有如下的有益效果：
[0087] 1)星敏感器是观测
恒星的天体敏感器，而利用星敏感器进行星间相对测量需要满足特定条件，本发明提出星间观测需要的光照条件和星敏感器观测条件，解决传统星敏感器只能被动观测问题，提高自主选星准确性；
[0088] 2)在实现星间观测
基础上，本发明提出实时计算主星相对子星方位矢量和方位角和俯仰角方法，并且利用万向轴调整星敏感器光轴指向连续
跟踪子星，解决传统观测无法连续跟踪问题，提高星间连续观测效率。
附图说明
[0089] 图1为本发明具体实施方式的方法
流程图；
[0090] 图2为本发明具体实施方式中主星星敏感器观测子星的流程图；
[0091] 图3为本发明具体实施方式中主星相对子星星间特定距离范围示意图；
[0092] 图4为本发明具体实施方式
中子星光照条件示意图；
[0093] 图5为本发明具体实施方式中星敏感器视场与地球位置关系示意图；
[0094] 图6为本发明具体实施方式
中子星可视星等计算示意图；
[0095] 图7为本发明具体实施方式中子星在星敏感器二维像面阵投影示意图；
[0096] 图8为本发明具体实施方式中主星相对子星方向矢量与方位角示意图。
具体实施方式
[0097] 下面结合具体实施方式和附图对本发明的技术方案作进一步的介绍。
[0098] 本具体实施方式公开了一种基于星敏感器和星间链路的HEO卫星编队飞行自主导航方法，在HEO卫星对地观测任务阶段，HEO卫星利用星敏感器和星间链路进行星间观测，是一种非常适合于HEO卫星自主连续相对观测的方法。如图1和图2所示，包括以下步骤：
[0099] S1：以HEO卫星对地观测为任务需求，将两个HEO卫星分别设为主星和子星，设计主星和子星编队飞行构型及轨道参数；轨道参数包括轨道半长轴a、轨道偏心率e、轨道倾角i、升交点赤经Ω、近地点幅角ω、过近地点时刻tp。
[0100] S2：根据地心惯性坐标系下卫星相对轨道动力学模型，建立自主导航系统状态模型。
[0101] S3：根据步骤S1计算主星和子星相对距离，如图3所示，判断子星是否满足星敏感器观测距离要求：如果满足，则进入步骤S4；否则，进入步骤S12。
[0102] S4：当主星观测子星时，子星需要被太阳光充分照射。当子星在地球光照区时，子星能被太阳光充分照射；反之，当子星进入地球阴影区时，由于地球遮挡，太阳光无法照射到子星，因此需要对子星光照条件进行判断。根据步骤S1解算太阳、地球和子星三者的位置关系，判断子星是否处在太阳光照区：如果是，则进入步骤S5；否则，进入步骤S12。
[0103] S5：根据步骤S1解算地球、主星和子星三者的位置关系，判断地球否进入星敏感器视场：如果是，则进入步骤S6；否则，进入步骤S12。
[0104] S6：根据步骤S1计算子星可视星等，判断子星可视星等是否小于星敏感器可观测阈值：如果是，则进入步骤S7；否则，进入步骤S12。
[0105] S7：根据步骤S1计算主星相对子星方向矢量与星敏感器光轴指向夹角，判断子星是否在星敏感器视场范围内：如果是，则进入步骤S8；否则，通过万向轴调整星敏感器光轴指向，再次按照步骤S7判断子星是否在星敏感器视场范围内，是则进入步骤S8，否则进入S12；
[0106] S8：根据步骤S7计算子星在星敏感器二维像面阵坐标，判断子星是否在星敏感器二维像面阵内：如果是，则进入步骤S9；否则，进入步骤S12。
[0107] S9：计算主星相对子星的理论方向矢量、方位角以及俯仰角。
[0108] S10：调整星敏感器真实光轴指向与理论方向矢量一致，对子星进行真实观测，计算子星相对卫星真实方向矢量，建立以单位方向矢量和距离为观测量的观测方程。
[0109] S11：对步骤S2所建立的状态模型和步骤S10所建立的观测方程离散化，利用Unscented卡尔曼滤波算法估计卫星位置和速度。
[0110] S12：结束观测。
[0111] 步骤S2中，建立自主导航系统状态模型的过程如下：
[0112] 在地心惯性坐标系下，当主星位置距离大于子星与主星相对距离的时候，建立卫星相对目标子星轨道动力学模型：
[0113]
[0114]
[0115] 其中，δr(10)为主星相对子星方向矢量，δv(10)为主星相对子星速度矢量， 为卫星位置随时间的一阶导数， 为卫星速度随时间的一阶导数，r(0)为主星位置矢量，μe为地球引力常数，af为摄动加速度；
[0116] 定义状态变量x＝[(δr(10))T(δv(10))T]T，建立自主导航系统状态模型；
[0117]
[0118] 其中，f为系统非线性连续状态转移函数，wt为状态噪声。
[0119] 步骤S3中，判断子星否满足星敏感器观测距离要求的过程如下：计算主星相对子星距离δr(10)，如图3所示，判断δr(10)是否满足式(3)所示的条件：
[0120] Lmin≤δr(10)≤Lmax (3)
[0121] 其中，δr(10)＝|δr(10)|＝|r(1)-r(0)|，r(0)为主星位置矢量，r(1)为子星位置矢量，Lmin(10) (10) (10)星间观测所需最小距离，Lmax为星间观测所需最大距离。δr 和δr 的区别在于：δr 用正常字体r来表示标量，而δr(10)用加粗字体来表示矢量r，这是本领域的常用表示方式，其他参数也同理，不再赘述。
[0122] 根据太阳、地球和子星三者几何位置关系，如图4所示，确定太阳阴影区和子星运行轨迹穿过该阴影区的临界条件。步骤S4中，判断子星是否处在太阳光照区的过程如下：分析地球阴影范围以及子星运行穿过该阴影范围的临界条件，设子星位置矢量r(1)与太阳位置矢量r(sun)夹角为ψ，子星进入和离开地球阴影范围的临界夹角为ψcri，则子星处在太阳光照区需要满足条件：
[0123] ψ＜ψcri (4)。
[0124] 其中， Re是地球半径。
[0125] 步骤S5中，判断地球是否进入星敏感器视场的过程如下：如图5所示，设主星位置矢量r(0)和主星相对子星方向矢量δr(10)的夹角为θ，被地球遮挡导致背景光线过弱的临界条件是主星相对子星方向矢量δr(10)与地球边缘相切，定义此临界夹角为θcri，则地球未进入星敏感器视场条件为：
[0126] θ＞θcri (5)。
[0127] 步骤S6中，子星可视星等根据式(6)计算得出：
[0128]
[0129] 式(6)中，m为子星可视星等，M为子星绝对星等，根据式(7)计算得到；|r(sun0)|是太(10)阳与子星之间的距离；如图6所示，ξ是主星相对子星方向矢量δr 与太阳相对子星方向矢量r(sun1)的夹角，p(ξ)是相位积分，根据式(8)计算得到；d0是地球与太阳之间的平均距离；
[0130]
[0131] 式(7)中，msun是太阳的可视星等，rd为被观测天体的半径，a是被观测天体的反射率；
[0132]
[0133] 步骤S7中，判断子星是否在星敏感器视场范围内的过程如下：如果 则子星在星敏感器视场范围内；否则，则子星不在星敏感器视场范围内；如式(9)所示，FOV为星敏感器视场角；
[0134]
[0135] 式(9)中，δr(10)为主星相对子星方向矢量， 为星敏感器方向矢量， 是地心惯性坐标系相对本体坐标系姿态转换矩阵；
[0136] 步骤S8中，判断子星是否在星敏感器二维像面阵内的过程如下：如果满足式(10a)和(10b)，则子星在星敏感器二维像面阵内；
[0137]
[0138]
[0139] 其中， 为主星相对子星方向矢量δr(10)投影在星敏感器二维像面阵的坐标，如图7所示，IPlongth为二维像面阵的长度，IPwidth为二维像面阵的宽度。
[0140] 步骤S9中，主星相对子星的理论方向矢量、方位角和俯仰角分别根据式(11)、(12a)和(12b)计算得到：
[0141]
[0142] 式(11)中，δr(10)为主星相对子星的理论方向矢量，也称为主星相对子星方向矢量； 为主星相对子星单位方向矢量，主星相对子星理论矢量方向可由方位角和俯仰角描述，如图8所示，在卫星本体坐标系ob-xbybzb中，定义方位角α为δr(10)在ob-ybzb平面的投影与yb轴夹角，俯仰角δ为δr(10)与xb轴夹角，可表示为；
[0143]
[0144]
[0145] 其中， 是地心惯性坐标系相对本体坐标系姿态转换矩阵，α为方位角，δ为俯仰角。
[0146] 步骤S10中，建立的观测方程为：
[0147]
[0148] 式(13)中， 为主星相对子星单位方向矢量的真实测量值， 为主星相对子星距离的真实测量值。
[0149] 步骤S11的具体过程如下：
[0150] 根据式(14a)对步骤S2所建立的状态模型进行离散化：
[0151]
[0152] 式(14a)中，xk∈RL为第k个状态变量，xk+1为第k+1个状态变量，wk∈N(0,Qk)为第k个过程噪声，Qk为系统噪声强度，f为系统非线性连续状态转移函数；
[0153] 根据式(14b)对步骤S10所建立的观测方程进行离散化：
[0154] yk＝g(xk)+vk (14b)
[0155] 式(14b)中，yk∈RM为第k个输出向量，g为观测方程，vk∈N(0,Rk)为第k个测量噪声，Rk为观测噪声强度，wk和vk不相关；
[0156] 利用Unscented卡尔曼滤波算法估计卫星位置和速度，步骤如下：
[0157] S11.1：根据式(15a)、(15b)和(15c)进行Unscented变换：
[0158]
[0159]
[0160]
[0161] 其中， 为xk的均值， 为xk的方差，λ＝α2(1+κ)-1为标量参数，α决定了sigma点在 周围的分布，κ为标量参数，L为sigma点维数，χi,k-1是中间变量；
[0162] S11.2：预测过程：根据式(16a)、(16b)、(16c)、(16d)、(16e)和(16f)进行预测；
[0163] χi,k/k-1＝f(χi,k-1) (16a)
[0164]
[0165]
[0166]
[0167]
[0168]
[0169] 其中，Wi和 分别为U变换计算均值和方差所用的权系数，χi,k/k-1、Yi,k/k-1和 都是中间变量；
[0170] S11.3：更新过程：根据式(17a)、(17b)、(17c)、(17d)和(17e)进行更新；
[0171]
[0172]
[0173]
[0174]
[0175]
[0176] 其中， 和Kk都是中间变量；
[0177] S11.4：返回步骤S11.1进行下一个周期的滤波。
[0178] 综上，本具体实施方式首先以HEO卫星对地观测为任务需求，设计两个卫星(设为主星和子星)编队飞行构型和轨道参数，然后根据地心惯性坐标系下主星相对子星轨道动力学模型，建立自主导航系统状态模型；其次提出主星星敏感器观测子星所需满足的理论光照条件和成像条件。计算主星相对子星理论方位角与俯仰角，调整主星星敏感器真实光轴与理论方向一致，对子星进行真实观测，建立以相对单位方向矢量和距离为观测量的观测方程；最后使用Unscented卡尔曼滤波估计卫星位置和速度，本具体实施方式能够有效修正卫星相对位置误差，提高相对导航精度，非常适用于卫星编队飞行自主导航，本具体实施方式属于航天导航技术领域，不仅可以为卫星编队飞行提供高精度导航信息，而且可以为其自主导航系统设计提供参考。
