技术领域
[0001] 本发明属于多源信息导航领域，具体而言，涉及一种基于非合作近地轨道常驻群体观测的自主导航方法。
相关背景技术
[0002] 在无人机导航技术方面，使用最为广泛和成熟的导航方式是GNSS/INS组合导航技术。但卫星导航存在容易受到电磁干扰的弊端，在未来战争高强度电子对抗模式下容易因为敌方干扰而失去作用，因此研究在卫星导航失效情况下可用的天文导航技术显得尤为迫切。天文导航是目前除惯性导航以外抗干扰能力最好的导航方式，这是因为天文导航的观测对象是无法毁灭的自然天体，具备安全可靠、被动探测不易暴露、无累积误差等战术优越性，因此受到各个国家的充分重视。近年来，随着人造卫星的增多，可观测的近地空间物体也越来越多，与只检测假定为“无穷远”的恒星不同，近地空间物体相对于观测者上传感器的位置随着轨道位置的变化而变化。这种轨道位置和姿态依赖性意味着，理论上，除了通常为恒星跟踪器假设的传统姿态信息外，恒星跟踪器还可以提供位置信息。因此，依托于混合识别恒星及近地空间物体的天文导航系统也开始发展。
[0003] 本发明提供了一种基于非合作近地轨道常驻群体观测的自主导航方法，依托于观测近地空间物体以获得载体的测量位置和姿态信息，该方法相比于传统的因子图融合方法，针对观测恒星定位定姿都受水平基准精准度影响的问题，以RSO观测传感器测量解算的载体位置和姿态构建量测方程，以及使用观测得到的RSO矢量和RSO位置信息构建量测方程。此外，该方法以惯性导航系统为核心传感器构建因子图，根据量测方程构建RSO传感器因子并加入因子图解算，以获得最大后验概率优化导航结果。该方法相比传统观测恒星的CNS/INS因子图导航信息融合方法，有效提升了位置优化的准确性。
具体实施方式
[0108] 下面结合附图对本发明的技术方案做进一步的详细说明，所述实施方式的示例在附图中示出，下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。
[0109] 本发明是一种基于非合作近地轨道常驻群体观测的松组合自主导航方法(方法1)，以及一种基于非合作近地轨道常驻群体观测的紧组合自主导航方法(方法2)，其方法流程分别如图1、图2所示，包括以下几个步骤：
[0110] 方法1、方法2的步骤1，采集惯性导航系统的角速度和加速度输出数据。所述步骤1的具体步骤包括以下步骤：
[0111] 步骤1.1，采集载体中惯性导航系统的加速度计测量数据，加速度计的输出数据如下公式：
[0112]n
[0113] 其中，f 为惯性导航系统的加速度计所输出的比力， 为载体相对于地球的加速度，ωie为地球自转角速度，ωen为载体相对于地球的角速度，Ven为载体相对于地球的运动速度，g为载体处的地球的重力加速度。
[0114] 步骤1.2，采集载体中惯性导航系统的陀螺仪测量数据，陀螺仪的输出数据如公式所示：
[0115]
[0116] 其中， 为载体系的陀螺仪相对于惯性系的角速度在载体系上的投影， 为地理系相对于惯性系的角速度在载体系上的投影， 为导航系相对于地理系的角速度在载体系上的投影， 为载体系相对于地理系的角速度在载体系上的投影；
[0117] 方法1、方法2的步骤2，采集观测RSOs在观测坐标系下的单位矢量以及RSOs位置信息，所述步骤2的具体步骤包括以下步骤：
[0118] 步骤2.1，星敏感器图像传感器拍摄当前视场范围内的RSO所成图像，对图像进行预处理，选择信噪比(SNR)大于用户定义的检测阈值的像素，并将其他像素丢弃。
[0119] 步骤2.2，由预处理后的图像，利用最小距离滤波器对两个聚类进行判断，最小距离滤波判断如公式所示：
[0120] δmin≤εdist
[0121] 其中δmin表示用一致范数定义的两个聚类之间的距离，εdist为用户定义的阈值，两个聚类满足最小距离滤波判断则通过最小距离滤波器，进行后续的判断。
[0122] 步骤2.3，如果两个聚类满足最小距离滤波判断，则利用增长长度滤波器进行判断，聚类Mi的长度向量li＝[ly,i,lz,i]表达式为：
[0123]
[0124] 定义逻辑变量L1、L2和L3，表达式为：
[0125]
[0126] 三个逻辑变量L1、L2和L3中至少有一个为1则增长长度滤波器成功通过，进行后续的判断。
[0127] 步骤2.4，如果两个聚类通过最小距离滤波器和增长长度滤波器，则利用密度滤波器进行判断，聚类Mi的密度定义如公式所示：
[0128]
[0129] 其中Ni是属于第i个聚类的像素数量，且λi的表达式为：
[0130]
[0131] λi表示一个虚拟长度，它足以定义密度而不需要复杂的计算。角度αi的表达式为：
[0132]
[0133] 两个聚类的密度相对差异表达式为：
[0134]
[0135] 其中di和dj是单个原始聚类的密度，dij是合并聚类的密度。两个聚类的密度相对差异必须低于用户定义的阈值εdens,rel，逻辑变量D1表达式为：
[0136]
[0137] 逻辑变量D2表达式为：
[0138]
[0139] 其中εdens,abs是用户定义的阈值。如果D1＝D2＝1，则密度过滤器成功通过。
[0140] 步骤2.5，如果两个聚类通过最小距离滤波器、增长长度滤波器和密度滤波器，则T可以进行合并，并计算合并后的质心。质心ci＝[cy,i,cz,i]是第i个聚类的参考位置，它是通T
过整个聚类Mi中的像素坐标pi＝[py,pz]上进行能量加权平均来计算的，表达式为：
[0141]
[0142] 其中E(pj)是像素pj的信号强度，第i个聚类的总信号强度表达式为：
[0143]
[0144] 步骤2.6，由星图识别算法在RSO星库中找到各聚类的对应匹配，获取各RSO星历信息。
[0145] 步骤2.7，将星敏感器坐标系与机体坐标系对齐，通过RSO在焦平面上成像的坐标RSO(xc,yc)，获得RSO矢量在星敏感器坐标系中表示的单位矢量s ，表达式如公式所示：
[0146]
[0147] 式中， 为RSO矢量与星敏感器坐标系所成的方位角和高度角；f为星敏感器光学透镜焦距。
[0148] 步骤2.8，根据RSO在大地坐标系下的星历信息，获得在大地坐标系下表示的RSO位置 其表达式如公式所示：
[0149]
[0150] 其中 分别表示RSO在大地坐标系下的xyz三轴位置。
[0151] 方法1的步骤3，根据测量的RSOs方位角、高度角以及位置信息解算出载体位置和姿态信息并构建量测方程，所述步骤3的具体步骤包括以下步骤：
[0152] 步骤3.1，设置导航坐标系n为原点在 的东北天地理坐标系，其对T
应的纬经度为[L,λ] ,根据RSO在大地坐标系下表示的位置 获得RSO在导航坐标系下表示的位置 其表达式如公式所示：
[0153]
[0154] 其中， 表示地球系e到导航系n的旋转矩阵。
[0155] 步骤3.2，计算星敏感器观测到的RSO单位矢量vRSO，其表达式如公式所示：
[0156]
[0157] 其中，pn＝[xn,yn,zn]T表示载体在导航坐标系中的位置。
[0158] 步骤3.3，设载体相对于导航系的姿态角为 计算导航系n到机体系b的旋转矩阵 其表达式如公式所示：
[0159]
[0160] 步骤3.4，根据RSO矢量在星敏感器坐标系下和在导航系下表示的单位矢量s和v,RSO RSO求得s 和v 之间的转换关系，其表达式如公式所示：
[0161]
[0162] 步骤3.5，对于观测到的第m个RSO，令观测到的RSO矢量为sRSOm，RSO在导航系下的位置为 则有：
[0163]RSOm
[0164] 其中s 、 已知，则通过多个观测到的RSO，可以解出 xn。
[0165] 步骤3.6，得到观测RSO解算的量测量 量测方程表达式如公式所示：
[0166]
[0167] 其中xn,yn,zn,γn,θn,ψn表示载体在地理系下的东北天位置、横滚角、俯仰角和航RSO向角真值，n 表示量测噪声。
[0168] 方法2的步骤3，直接根据测量的RSOs方位角、高度角以及位置信息推算出RSO量测方程，所述步骤3的具体步骤包括以下步骤：
[0169] 步骤3.1，设置导航坐标系n为原点在 的东北天地理坐标系，其对T
应的纬经度为[L,λ] ,根据RSO在大地坐标系下表示的位置 获得RSO在导航坐标系下表示的位置 其表达式如公式所示：
[0170]
[0171] 其中， 表示地球系e到导航系n的旋转矩阵。
[0172] 步骤3.2，计算星敏感器观测到的RSO单位矢量vRSO，其表达式如公式所示：
[0173]
[0174] 其中，pn＝[xn,yn,zn]T表示载体在导航坐标系中的位置。
[0175] 步骤3.3，设载体相对于导航系的姿态角为 计算导航系n到机体系b的旋转矩阵 其表达式如公式所示：
[0176]
[0177] 步骤3.4，根据RSO矢量在星敏感器坐标系下和在导航系下表示的单位矢量s和v,RSO RSO求得s 和v 之间的转换关系，其表达式如公式所示：
[0178]
[0179] 步骤3.5，设sRSO为观测RSO的量测值，计算sRSO对状态量[xn,yn,zn,γn,θn,ψn]T的部分观测矩阵H，其表达式如公式所示：
[0180]
[0181] 步骤3.6，得到观测RSO的量测方程，表达式如公式所示：
[0182]
[0183] 其中nsRSO表示量测噪声。
[0184] 方法1、方法2的步骤4，由惯性导航系统与RSO观测传感器构建因子图，所述步骤4的具体步骤包括以下步骤：
[0185] 步骤4.1，利用步骤2所提供的RSO观测传感器测量数据，根据步骤3所述过程构建RSO传感器因子，其量测方程 简写如公式所示：
[0186]
[0187] 其中，hRSO(·)是tk时刻RSO观测传感器的量测模型，nRSO为RSO观测传感器的测量噪声。则RSO传感器因子如公式所示：
[0188]
[0189] 其中，d(·)为代价函数，err(·)为误差函数。
[0190] 步骤4.2，在有新的RSO观测传感器因子建立后，根据步骤1采集到的惯性导航系统的角速度和加速度测量值，构建IMU因子，为了避免在IMU数据过多时建立因子，造成系统计算压力过大，采用预积分形式的IMU因子，其思想是将因子相关的两个时间实例tk和tk+1之间的连续IMU测量集成到单个分量中。由此得到IMU预积分因子如公式所示：
[0191]
[0192] 其中，xk+1和xk分别表示tk+1时刻与tk时刻的位置、速度与姿态状态量，αk表示惯性导航器件在tk时刻的偏差变量，Δxk→k+1表示包括从tk时刻到tk+1时刻载体位置、速度和姿态的累积变化，h(xk,αk,zk)表示由状态方程得到的tk+1时刻的导航状态量预测值，其与当前估计值xk+1的差值即为IMU预积分因子代表的需要进行的最小化误差函数，d(·)为代价函数。
[0193] 方法1、方法2的步骤5，利用因子图模型对代价函数进行分解，以使得每一个因子节点都跟传感器对应，由此得到简化后的导航状态最大后验估计如公式：
[0194]
[0195] 其中，fi(·)为因子节点对应的局部函数，Xi为对应传感器因子节点的状态量集合，X为导航系统所有状态变量节点的集合；
[0196] 导航状态量和各传感量测值的后验概率密度函数p(X|Z)与所有因子节点乘积成正比，即有公式：
[0197]
[0198] 其中，X是状态变量的集合，Z是传感器量测值集合， 为马氏距离的平方，Σi为对应的协方差矩阵，zi是因子节点θi处的量测值。由此，各传感器的融合问题就是求解最大后验估计 即为方法所求的最优估计，而系统状态量的最大后验估计本质上就是求解下列非线性最小二乘公式，如公式：
[0199]
[0200] 实施例：
[0201] 为了验证本发明所提出的方法1，即一种基于非合作近地轨道常驻群体观测的自主导航方法的有效性，进行数字仿真分析。仿真以惯导为核心传感器，RSO观测传感器为辅助传感器，使用仿真数据得到的惯导数据与RSO观测传感器测角数据进行组合导航算法的验证。
[0202] 使用本方法进行RSO观测惯性系姿态角观测下的因子图优化结果作为实例说明。使用仿真数据的验证结果如图3，4，5，6，7，8，9所示。仿真数据采用一条设定的2000秒机载轨迹，载体观测的RSOs星历采用STARLINK星链星历信息，观测时间为UTC时间2024年05月28日02点20分00秒开始的2000s的数据，并设定星敏感器观测RSO的方位角和高度角误差为20角秒；本发明方法1和真值数据作差的位置误差曲线如图4所示，本发明方法1和真值数据作差的速度误差曲线如图5所示，本发明方法1和真值数据作差的姿态角误差曲线如图6所示，本发明方法1、传统惯性/天文导航方法和真值数据作差的位置误差对比如图7所示，本发明方法1、传统方法和真值数据作差的速度误差对比如图8所示，本发明方法1、传统方法和真值数据作差的姿态角误差对比如图9所示。
[0203] 由图4，5，6可以看出，本发明方法1通过观测RSO解算可以对位置速度姿态同时进行优化，优化后位置误差都保持在20m内，速度误差基本保持在0.4m/s内，最高时0.4m/s，姿态误差基本保持在20角秒之内。
[0204] 由图7，8，9可以看出，传统惯性/天文导航方法主要对姿态角有修正作用，优化后姿态误差基本保持在20角秒之内，对于位置和速度的修正效果则很小，而与其相比，本发明方法2对于位置速度姿态的优化效果均有极大提升，解决了传统天文导航无法单独修正位置的问题。
[0205] 为了验证本发明所提出的方法2，即一种基于观测近地空间物体的因子图导航信息融合方法的有效性，进行数字仿真分析。仿真以惯导为核心传感器，RSO观测传感器为辅助传感器，使用仿真数据得到的惯导数据与RSO观测传感器测角数据进行组合导航算法的验证。
[0206] 使用本发明方法2进行观测RSOs的因子图优化结果作为实例说明。使用仿真数据的验证结果如图3，10，11，12，13，14，15，16，17，18，19所示。仿真数据与发明方法1同样采用图3设定的2000秒机载轨迹，载体观测的RSOs星历采用星链星历信息，观测时间为UTC时间2024年05月28日02点20分00秒开始的2000s的数据，并设定星敏感器观测RSO的方位角和高度角误差为20角秒；本发明方法2和真值数据作差的位置误差曲线如图10所示，本发明方法
2和真值数据作差的速度误差曲线如图11所示，本发明方法2和真值数据作差的姿态角误差曲线如图12所示，本发明方法2、传统惯性/天文导航方法和真值数据作差的位置误差对比如图13所示，本发明方法2、传统方法和真值数据作差的速度误差对比如图14所示，本发明方法2、传统方法和真值数据作差的姿态角误差对比如图15所示。仿真时间内RSOs的运行轨迹如图16所示，仿真时间内在导航坐标系内可观测到的RSOs运行轨迹如图17所示，可观测的RSOs数量随时间的变化如图18所示，仿真时间内星敏感器坐标下STARLINK‑31359的方位角和高度角随时间变化如图19所示。
[0207] 由图10，11，12可以看出，本发明方法2通过观测RSO解算可以对位置速度姿态同时进行优化，优化后位置误差都保持在13m内，速度误差基本保持在0.2m/s内，最高时0.4m/s，姿态误差基本保持在10角秒之内。
[0208] 由图13，14，15可以看出，传统惯性/天文导航方法主要对姿态角有修正作用，优化后姿态误差基本保持在20角秒之内，对于位置和速度的修正效果则很小，而与其相比，本发明方法2对于位置速度姿态的优化效果均有极大提升，解决了传统天文导航无法单独修正位置的问题。
[0209] 本技术领域技术人员可以理解，除非特意声明，这里使用的单数形式“一”、“一个”、“所述”和“该”也可包括复数形式。应该进一步理解的是，本发明的说明书中使用的措辞“包括”是指存在所述特征、整数、步骤、操作、元件和/或组件，但是并不排除存在或添加一个或多个其他特征、整数、步骤、操作、元件、组件和/或它们的组。应该理解，当我们称元件被“连接”或“耦接”到另一元件时，它可以直接连接或耦接到其他元件，或者也可以存在中间元件。此外，这里使用的“连接”或“耦接”可以包括无线连接或耦接。这里使用的措辞“和/或”包括一个或更多个相关联的列出项的任一单元和全部组合。
[0210] 本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。
[0211] 本技术领域技术人员可以理解的是，本发明可以涉及用于执行本申请中所述操作中的一项或多项操作的设备。所述设备可以为所需的目的而专门设计和制造，或者也可以包括通用计算机中的已知设备，所述通用计算机有存储在其内的程序选择性地激活或重构。这样的计算机程序可以被存储在设备(例如，计算机)可读介质中或者存储在适于存储电子指令并分别耦联到总线的任何类型的介质中，所述计算机可读介质包括但不限于任何类型的盘(包括软盘、硬盘、光盘、CD‑ROM、和磁光盘)、随机存储器(RAM)、只读存储器(ROM)、电可编程ROM、电可擦ROM(EPROM)、电可擦除可编程ROM(EEPROM)、闪存、磁性卡片或光线卡片。可读介质包括用于以由设备(例如，计算机)可读的形式存储或传输信息的任何机构。例如，可读介质包括随机存储器(RAM)、只读存储器(ROM)、磁盘存储介质、光学存储介质、闪存装置、以电的、光的、声的或其他的形式传播的信号(例如载波、红外信号、数字信号)等。
[0212] 本技术领域技术人员可以理解的是，可以用计算机程序指令来实现这些结构图和/或框图和/或流图中的每个框以及这些结构图和/或框图和/或流图中的框的组合。可以将这些计算机程序指令提供给通用计算机、专业计算机或其他可编程数据处理方法的处理器来生成机器，从而通过计算机或其他可编程数据处理方法的处理器来执行的指令创建了用于实现结构图和/或框图和/或流图的框或多个框中指定的方法。
[0213] 本技术领域技术人员可以理解的是，本发明中已经讨论过的各种操作、方法、流程中的步骤、措施、方案可以被交替、更改、组合或删除。进一步地，具有本发明中已经讨论过的各种操作、方法、流程中的其他步骤、措施、方案也可以被交替、更改、重排、分解、组合或删除。进一步地，现有技术中的具有与本发明中公开的各种操作、方法、流程中的步骤、措施、方案也可以被交替、更改、重排、分解、组合或删除。
[0214] 以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。
