CN112161632A - 一种基于相对位置矢量测量的卫星编队初始定位算法 - Google Patents

Abstract

本发明提出了一种基于相对位置矢量测量的卫星编队初始定位算法。本算法通过将星间相对加速度与引力场模型相结合，反解得到卫星编队的绝对位置。开发了通过拟合相对测量数据，得到星间相对加速度的自适应多项式拟合算法。给出了不同导航解算模式选取的判别依据。导航结果可作为后续精密定轨的卫星状态量初值。本发明完全不依赖外部信息，具有自主性高、适应能力强的优点。

Description

一种基于相对位置矢量测量的卫星编队初始定位算法

技术领域

本发明提供一种具有完全自主性的、不需要外部测量信息的卫星编队初始定位算法。它涉及一种从一系列星间相对位置矢量中提取相对加速度的自适应拟合阶次的多项式拟合算法，以及基于星间相对运动特性的导航解算算法，包括导航解算模式选取判据和全局/局部搜索算法。属于导航技术领域。

背景技术

卫星编队(也称为分布式卫星系统)作为航天领域的一个创新概念，代表了航天技术研究的一个重要方向。卫星编队与单颗卫星相比，由于其卫星数量更多，在功能和性能上得到了提高，并且具有极大的应用灵活性和安全性，成为当下研究的热点。

随着各国相关卫星编队发射任务的增加，在轨运行卫星数量的增加对传统依靠地面站测控获得卫星精密轨道数据的方法带来越来越大的挑战，因此对卫星自主导航方法的需求在日益增长。目前低轨卫星主要利用全球导航卫星系统 GNSS(Global NavigationSatellite System)完成自主导航。但是如果在战时地面站和GNSS卫星全部被摧毁或者是失效的情况下，依靠地面站和GNSS的导航方法都将失效。并且随着深空探测任务的发展，深空探测器需要在没有 GNSS和地面站信息的情况下完成自主导航。因此研究一种完全独立的卫星编队自主导航方法具有重要意义。

作为一种新型自主导航方法，基于星间相对测量的卫星编队自主导航方法可以满足上述需求。

目前，很多研究已经证明了利用相对位置测量实现两个或两个以上编队航天器自主定轨的可行性。这种自主导航方法的基本原理是通过卫星编队卫星对彼此之间的相对位置矢量进行一段时间的连续测量，进而对所有成员航天器进行绝对轨道确定，其中星间相对距离可以用微波测距系统或激光测距仪测量，星间相对视线(LOS)由光学照相机来测量。此外为了确定卫星编队内各卫星的绝对轨道，各卫星还需要配备星敏感器以确定卫星在惯性空间中的姿态从而将星间相对位置矢量转换到惯性空间内。

基于星间相对位置测量确定绝对轨道的自主导航方法最早由Muller和 Kachmar在1971年提出。他们在论文中提出了一种“卫星弹射”导航系统，即主卫星向临近空间中弹射一颗小卫星。并且主卫星可以利用星上设备跟踪测量小卫星相对于主卫星的位置变化。一段时间内的连续相对位置矢量测量信息经过序贯卡尔曼滤波器处理，进而可以得到对主卫星绝对轨道状态量的估计值，实现自主导航。

1983年，K.Yong,C.C.Chao和A.S.Liu提出利用星间激光通信进行星间相对位置矢量测量，进而确定卫星编队绝对轨道的方法。该方法首先利用卡尔曼滤波器对相对测量进行估计得到精确测量信息，之后利用批处理滤波器对精确相对测量信息进行处理，确定卫星编队绝对轨道。仿真结果表明，定位结果轨道面内误差为200m，轨道面外误差为300m。

Markley在1984年通过将系统模型进行线性化对这种导航方法的可观测性进了分析。他指出几乎所有的双卫星编队都可以通过相对位置矢量测量数据确定绝对轨道，除非两颗卫星的轨道高度始终相同。

1999年，Psiaki给出了这种自主导航方法原理的一种解释。他指出卫星编队可以被看作为一个长基线重力梯度探测仪，其中每一颗卫星为加速度计的质量块。星间相对位置矢量变化是由卫星编队中心点处的重力梯度引起的，这使得可以确定整个卫星编队的绝对轨道。Psiaki在之后研究中提出了一种可以在确定卫星绝对轨道的同时修正中心天体重力场模型的算法。他在论文中对一个绕月双星编队进行了仿真，仿真结果表明该算法得到的绝对位置均方根误差为 10m，绝对速度均方根误差为0.01m/s，以及重力场模型均方根误差为2mGal(2 ×10^-5m/s²)。

最新研究结果也证明了基于星间相对测量的自主导航方法的有效性。ChrisR.Gnam，Andrew D.Dianetti和John L.Crassidis利用星间相对视线(LOS)和星间距离测量对多星编队的轨道确定问题进行了分析。仿真结果表明，通过使用UKF滤波器和URTS平滑对相对测量数据进行处理，对于各种不同的仿真场景，定位结果均能保持在100m以内。

除了利用完整相对位置信息，即相对距离和相对视线构成的相对位置矢量信息，确定卫星集群绝对轨道以外，还有关于只利用部分相对位置信息进行绝对轨道确定的研究。Jo Ryeong Yim，John L.Crassidis和John L.Junkins研究了仅利用星间相对视线(LOS)测量对双星编队进行绝对定轨的可行性。仿真结果表明定位精度为200m，定速精度为0.2m/s。

1998年，P.A.M.Abusali,B.D.Tapley和B.E.Schutz指出GPS BlockⅡR 卫星可以进行星间测距。通过星间测距，GPS卫星可以向用户提供精度更好的卫星状态信息。仿真结果表明，通过星间距离测量可以有效抑制由GPS卫星长期轨道预报动力学模型和初值引起的误差，但是不能修正由地球自转参数(EOP)预报引起的误差。

2007年，Keric Hill和George H.Born证明了仅利用星间相对测距对运行于地月或地日L1、L2拉格朗日点HALO轨道上的卫星编队进行自主导航的可行性。仿真结果表明，导航误差在数十米量级。其后续研究表明在二体问题情况下利用相对距离信息只能确定相对轨道，但是当三体引力占主导时，仅利用星间距离测量就可以确定出卫星编队的绝对和相对轨道。利用EKF对测量数据进行处理，得到仿真结果为:对于运行于月球平动点附近轨道的卫星编队，定位均方根误差为量级为100m；对于运行于月球轨道的卫星编队，定位均方根误差量级为10m。

2012年，Jason M.Leonard和Brandon A.Jonesy等人通过对一个围绕小行星433Eros运行的双星编队进行仿真，得出即使中心天体的重力场模型未知，双星编队也可以通过星间相对距离和相对距离变化率测量实现绝对定轨。文中采用了一种迭代扩展卡尔曼平方根信息滤波器(IEK-SRIF)对相对测量信息进行处理，得到的仿真结果为，在长达9天的观测过程中，卫星编队位置均方根误差为56m，速度均方根误差为0.006m/s。

目前相关研究所采用的处理方法一般是利用卡尔曼滤波器或是批处理滤波器对星间相对测量数据进行处理，进而得到较高精度的卫星编队定轨/定位结果。而这些滤波方法，都需要一个相对准确的系统状态量初值。作为一种力图实现完全独立的、不需要任何外部信息的自主导航方法，基于星间相对测量的导航方法应该具备在没有任何先验信息的情况下，依然可以完成自主导航任务的能力。本发明提供了一种具有完全自主性的、不需要外部测量信息的卫星编队初始定位算法。该算法可以为后续基于相对测量的精密定轨算法提供所需的系统状态量初值。

发明内容

(一)发明目的

本发明为解决卫星编队初始定位问题，提出一种基于相对位置矢量测量的卫星编队初始定位算法。本发明具有完全自主性、不需要外部测量信息的优点，适用于卫星编队绝对导航定位的应用场景。

(二)技术方案

为了方便描述，以三颗星A、B、C构成的卫星编队为例，本发明所述的卫星编队初始定位算法实施步骤如下：

步骤一：通过自适应多项式拟合对相对位置变化进行表征

本发明采用对一系列相对位置矢量测量数据进行多项式拟合的方法表征星间相对运动，进而获取星间相对加速度信息。而作为一种完全自主的卫星编队导航算法，本发明需要在完全不依赖外部信息的情况下，根据测量数据的性质，自主选取对后续导航解算最有利的多项式拟合阶次，该过程称为自适应拟合阶次选取。

若拟合阶次较低，则出现欠拟合现象。此时模型拟合程度不高，数据距离拟合曲线较远。若拟合阶次过高，则会出现过拟合现象。此时拟合曲线受测量数据影响显著，偏离真实情况。因此要想更好地表征星间相对运动，获得精确的星间相对加速度信息，需要选取合适的拟合阶次。

本发明利用不同拟合阶次间拟合残差平方和的变化趋势，给出用于拟合阶次选取的自适应拟合判据。

首先，假设对一系列相对位置矢量测量数据进行n阶拟合，得到n阶相对运动拟合多项式为：

其中r_BA(t)表示t时刻卫星B与卫星A之间的相对位置矢量拟合值，下标x,y,z 分别表示相对位置矢量分量。α表示n阶多项式拟合的系数。

假设共有m个历元时刻的相对位置矢量测量数据，其测量模型可以表示为：

其中

和

分别代表星间相对距离和相对视线(LOS)的测量量。v_ρ代表星间相对距离的测量误差，认为其满足均值为零，标准差为σ_ρ的高斯分布。v_u为3×1的矢量，各分量分别代表相对视线的各轴测量误差，且认为三轴测量误差均满足均值为零，标准差为σ_u的高斯分布。

进而可得m次测量的测量协方差矩阵为：

需要指出的是，为了防止在拟合过程中出现数值问题，在拟合前首先需要对各观测历元时刻进行数值上的缩放：

进而通过对m个历元时刻的相对位置矢量测量进行加权最小二乘拟合，可得各阶拟合系数为：

P＝(H^TR^-1H))^-1

其中P为拟合系数的协方差矩阵，矩阵H为：

此外还可得到拟合残差平方和为：

本发明以相邻拟合阶次间残差平方和变化率作为拟合阶次选取指标，当残差平方和变化率小于一个小量ε时，此时阶次n即为最终自适应拟合算法选取的拟合阶次。具体定义方法如式(8)所示：

步骤二：星间相对加速度提取

根据步骤一选取得到的拟合阶次以及各阶拟合系数，得到星间相对运动拟合多项式如式(1)所示。求取式(1)的二阶导数，得到星间相对加速度多项式：

由此，可计算得到任一历元时刻的星间相对加速度矢量。为保证后续导航解算的精度，此处选取用于拟合的m个历元时刻的中间时刻作为导航解算时刻。利用式(9)计算得到导航解算时刻的星间相对加速度矢量，并从测量数据中提取出对应时刻的星间相对位置矢量测量数据。

步骤三：导航解算模式选取

卫星编队之所以可以通过连续的相对位置矢量测量完成卫星编队的绝对导航，是因为对于卫星编队来说，由于星间距离较近，所以如大气阻力、太阳光压和三体引力等摄动力可认为被抵消掉，而影响星间相对运动的只有中心天体引力。因此可通过将步骤二中得到的相对加速度和相对位置矢量与中心天体引力场模型相结合，反解得到卫星编队的绝对位置。

导航解算算法需要卫星编队位置的初始信息。在完全不依赖外部信息输入的情况下，卫星编队首先需要进行全局粗搜索，找到卫星编队位置的初始点。在这一过程中需要利用二体引力场模型。然而当两颗卫星始终具有相同轨道高度时，利用二体引力场模型解算得到的卫星编队位置会出现奇异，即满足相同相对位置矢量和相对加速度条件的点有无数个，无法完成全局粗搜索，此时需要新的导航解算模式。

导航解算模式选取的本质是要区分出不同卫星是否在导航解算时刻具有相同的轨道高度。本发明根据二体引力场模型的相关性质，给出利用相对位置矢量和相对加速度矢量之间的位置关系进行导航解算模式选取的判别依据。

当两颗卫星具有相同轨道高度时，星间相对位置矢量与相对加速度矢量反向，即：

因此给出导航解算模式选取判别依据为：

其中δ为一个大于-1的负数。当相对位置矢量与相对加速度矢量满足式(11)时，选取导航解算模式一，反之选取导航解算模式二。

步骤四：全局粗搜索算法

导航解算的本质是求解非线性方程组，但是典型的求解非线性方程组的方法都需要一个初始解，而通常初始解由外部输入或先验信息给出。本发明作为一种完全自主的导航方法，根据不同的导航解算模式给出了完全不需要外部输入和先验信息的、不同的初始解获取算法，由于获得的初始解只是对真实解在一定范围内的估计，故此部分算法称为全局粗搜索算法。

对于导航解算模式一，此时卫星编队的空间构型使得星间相对位置矢量与相对加速度矢量的夹角远离180°。上文已经指出，由于卫星编队星间距离较近，所以如大气阻力、太阳光压和三体引力等摄动力可认为被抵消掉，而影响星间相对运动的只有中心天体引力，即相对加速度主要由中心天体引力引起，对于二体引力模型，相对加速度如式(12)所示：

其中μ为地球引力常量，

分别表示星间相对位置矢量和相对加速度。

Psiaki的相关研究指出，星间相对运动由星间连线中点处的重力梯度引起，因此相对加速度可进一步简化为：

其中Γ(r_c)表示星间连线中点处的重力梯度。

将式(13)中的重力梯度矩阵转换为球坐标系形式为：

其中λ,θ分别代表中点的地心纬度和地心经度。通过将重力梯度矩阵写为球坐标形式，实现地心距与经纬度的解耦，进而可将问题由三维简化到二维。将求解式(13)转换为求目标函数极值问题。根据球坐标形式重力梯度矩阵，将式(13)转化为：

其中λ^*,θ^*分别为式(13)解的地心纬度和地心经度。

要解决如式(15)所示的求极值问题，首先采用二维二分法进行初步粗搜索。所谓二维二分法，即通过将经纬度球面划分为等间距的网格，寻找网格节点中的极小点。之后在极小点与相邻节点构成的局部区域内继续缩小节点间距，生成新的节点，进而搜索得到新的极小点。由此迭代若干步，可得到初步粗搜索结果。

初步粗搜索结果可作为后续优化算法的初始点。本发明采用Nelder-Mead 算法对目标函数极小点进行进一步搜索。该算法无需计算目标函数的导数信息，通过比较函数值对极值点进行搜索，是一种成熟的优化算法，具体算法流程此处不再赘述。在求得地心经纬度之后，地心距可由式(13-14)联立得到。实际上，式(13)有两个中心对称的解，其中错误的一个解可在后续轨道确定过程中剔除，并不影响本发明的使用。

对于导航解算模式二，卫星编队的空间构型使得星间相对位置矢量与相对加速度矢量的夹角接近180°。此时利用求极值的方法求解式(13)将出现数值困难，目标函数难以收敛到极值点附近，甚至当相对位置矢量与相对加速度矢量共线时，式(13)有无数解。因此需要采用新的全局粗搜索算法。

事实上，当相对位置矢量与相对加速度矢量共线时，式(13)的解是一个垂直于相对位置矢量的圆。由这一性质可进一步得到，对于同一卫星，若存在两个这样的、法向量不共线的圆，卫星的位置即是两圆的一个交点，因此对于卫星轨道高度相同的卫星编队，至少需要三颗星间相对位置矢量不共线的卫星才可以完成导航解算。此时卫星的位置矢量可通过两圆法向量的叉乘得到，如式(16)所示：

其中

为卫星A位置矢量的单位矢量的估计值，而位置矢量模长的估计值，即卫星A的地心距可由式(17)计算得到。

进而，卫星A位置的估计值为：

步骤五：局部精确搜索算法

全局粗搜索算法基于简化模型对卫星编队的位置进行初步估计，因此定位误差较大。要获得卫星编队的精确位置，还需要根据精确引力场模型构建高精度导航方程组。对于由三颗卫星A、B、C构成的卫星编队，导航方程组可以表述为：

其中

为星间相对位置矢量测量值，由星载敏感器测量得到，为已知量。

为星间相对加速度，由步骤二提取得到，为已知量。g为引力加速度，通过高阶引力场模型计算得到。

因此卫星编队的位置可通过利用牛顿迭代法求解导航方程组得到，牛顿迭代法相关计算公式如式(20-22)所示：

F(r_A)＝a-g (20)

其中

将全局粗搜索算法结果作为r_A的初始点

得到迭代公式为：

其中J为F关于r_A的偏导数矩阵：

其中Γ为各点处的重力梯度矩阵，由引力场模型对位置求偏导得到。

至此完成对卫星编队位置的导航解算。

(三)优点

本发明提供的基于相对位置矢量测量的卫星编队初始定位算法优点在于：

①本发明提出的算法具有完全自主性，可以在完全不需要外部信息的情况下实现对围绕中心天体运行的卫星编队的导航定位。适用于深空探测、探月等卫星与地面通信难度较大的应用场景。同时还适用于战时敌方对我方地面站和导航卫星进行干扰破坏，使得卫星编队难于接收到外部信息的应用场景。

②本发明提出的算法具有全面性，针对可能存在的不同卫星编队构型，给出了不同的导航解算算法，具有较强的适应能力。

附图说明

图1是本发明的实施步骤流程图

图2是全局/局部搜索算法流程图

具体实施方式

下面将结合附图1和技术方案，以三颗星A、B、C构成的卫星集群为例，对本发明的具体实施过程做进一步的详细说明。

步骤一：通过自适应多项式拟合对相对位置变化进行表征

本发明采用对一系列相对位置矢量测量数据进行多项式拟合的方法表征星间相对运动，进而获取星间相对加速度信息。通过分析不同拟合阶次间拟合残差平方和的变化趋势，本发明给出用于自适应拟合阶次选取的自适应拟合判据。其中拟合残差平方和定义式为：

根据相邻拟合阶次间残差平方和变化率作为拟合阶次选取指标，当残差平方和变化率小于一个小量ε时。此时阶次n即为最终自适应拟合算法选取的拟合阶次。自适应拟合判据的具体定义方法为：

本步骤对应附图1中的第一个方框。

步骤二：星间相对加速度提取

在得到各阶拟合系数之后，可以得到星间相对运动多项式。对该多项式求取二阶导数，得到星间相对加速度，进而实现星间相对加速度提取。为保证后续导航解算的精度，此处选取用于拟合的m个历元时刻的中间时刻作为导航解算时刻。

本步骤对应附图1中的第二个方框。

步骤三：导航解算模式选取

卫星编队构型不同，导航解算算法不同。本发明根据二体引力场模型的相关性质，给出了以相对位置矢量和相对加速度矢量之间的位置关系作为导航解算模式选取判别依据的判别方法。

当两颗卫星具有相同轨道高度时，星间相对位置矢量与相对加速度矢量反向，即：

因此给出导航解算模式选取判别依据为：

其中δ为一个大于-1的负数。当相对位置矢量与相对加速度矢量满足式(26)时，选取导航解算模式一，反之选取导航解算模式二。

本步骤对应附图1中的第三个方框。

步骤四：全局粗搜索算法

本发明作为一种完全自主的导航方法，根据不同的导航解算模式给出了完全不需要外部输入和先验信息的、不同的初始解获取算法。

对于导航解算模式一，此时卫星编队的空间构型使得星间相对位置矢量与相对加速度矢量的夹角远离180°。

首先根据Psiaki的相关研究，将相对加速度简化为：

其中Γ(r_c)表示星间连线中点处的重力梯度。将式(27)中的重力梯度矩阵转换为球坐标系形式为：

其中λ,θ分别代表中点的地心纬度和地心经度。通过将重力梯度矩阵写为球坐标形式，实现地心距与经纬度的解耦，进而可将问题由三维简化到二维。

之后，将求解式(27)转换为求目标函数极值问题：

其中λ^*,θ^*分别为式(27)解的地心纬度和地心经度。

要解决式(29)所示的求极值问题，首先采用二维二分法进行初步粗搜索。所谓二维二分法，即通过将经纬度球面划分为等间距的网格，寻找网格节点中的极小点。之后在极小点与相邻节点构成的局部区域内继续缩小节点间距，生成新的节点，进而搜索得到新的极小点。由此迭代若干步，可得到初步全局粗搜索结果。

之后以初步粗搜索结果作为初始点，采用Nelder-Mead算法对目标函数极小点进行进一步搜索。

最后，在求得地心经纬度之后，地心距可由式(27-28)联立得到。实际上，式(27)有两个中心对称的解，其中错误的一个解可在后续轨道确定中剔除。

对于导航解算模式二，卫星编队的空间构型使得星间相对位置矢量与相对加速度矢量的夹角接近180°。此时卫星的位置矢量可通过两个不共线的相对位置矢量叉乘得到，如式(30)所示：

其中

为卫星A位置矢量的单位矢量的估计值，而位置矢量模长的估计值，即卫星A的地心距可由式(31)计算得到。

进而，得到卫星A位置的估计值为：

本步骤对应附图1中的第四个方框，全局粗搜索算法计算流程图见附图2。

步骤五：局部精确搜索算法

首先根据中心天体引力场模型，构建导航方程组：

其中

为星间相对位置矢量测量值，由星载敏感器测量得到，为已知量。

为星间相对加速度，由步骤二提取得到，为已知量。g为引力加速度，通过高阶引力场模型计算得到。

之后利用牛顿迭代法求解导航方程组得到卫星编队位置，牛顿迭代法计算公式如式(34-36)所示：

F(r_A)＝a-g (34)

其中

将全局粗搜索算法结果作为r_A的初始点

得到迭代公式为：

其中J为F关于r_A的偏导数矩阵：

其中Γ为各点处的重力梯度矩阵，由引力场模型对位置求偏导得到。

至此完成对卫星编队位置的导航解算。本步骤对应附图1中的第五个方框，局部精确搜索算法计算流程图见附图2。

通过以上步骤，可以获得基于相对位置矢量测量的卫星编队初始定位结果。从而可以满足卫星编队在完全自主的情况下，仅利用相对位置矢量测量信息，实现绝对导航的需求。

Claims (1)

1.一种基于相对位置矢量测量的卫星编队初始定位算法，其特征在于：

步骤一：通过自适应多项式拟合对相对位置变化进行表征

通过将拟合残差平方和的变化趋势与自适应拟合判据比较，选取多项式拟合的最优阶次，进而得到可以最好地反映星间相对运动三个分量的拟合多项式。

步骤二：星间相对加速度提取

通过求取三个星间相对运动多项式的二阶导数，得到星间相对加速度矢量

步骤三：导航解算模式选取

卫星编队构型不同，全局/局部搜索算法也不同。因此在进行全局/局部搜索之前，需要完成对不同导航解算模式的选取。本发明给出了根据相对加速度矢量和相对位置矢量点积大小进行模式选取的导航模式选取判据。

步骤四：全局粗搜索算法

本发明包含两种导航解算模式，因此有两种完全不同的全局粗搜索算法。

当卫星编队中卫星的轨道高度不相等时，可由模式一进行导航解算。首先通过二维二分法得到初步全局粗搜索结果。之后利用Nelder-Mead算法搜索目标函数的极值点，得到最终全局搜索结果。

当卫星编队中卫星的轨道高度相等时，应选用模式二进行导航解算。此时卫星编队至少要包含三颗卫星，且星间相对位置矢量应满足不共线的要求。进而，全局粗搜索结果可由任意两个不共线的星间相对位置矢量叉乘得到。

步骤五：局部精确搜索算法

根据相对位置矢量、引力加速度和相对加速度的物理关系，构建导航方程组。并以步骤四中得到的全局粗搜索结果作为初始点，利用牛顿迭代法求解导航方程组。进而得到卫星编队的精确导航结果，该结果可作为后续精密定轨中卫星状态量的初始值。

Images