一种基于动态脉冲累积窗口的X射线脉冲星导航方法与流程

本发明涉及一种基于动态脉冲累积窗口的x射线脉冲星导航方法，属于航天器制导、导航与控制领域。

背景技术：

x射线脉冲星能够发射具有稳定周期的x射线脉冲信号，是理想的导航天体源，被誉为“宇宙中的灯塔”、“自然界的天体钟”。航天器上的脉冲到达时间观测量反映了航天器惯性位置矢量在脉冲星视线矢量上的投影，可用于确定航天器的位置和速度。近年来，x射线脉冲星导航已成为航天器自主导航领域的研究热点，理论方法不断完善。与此同时，x射线探测器技术也有了日新月异的发展，在探测灵敏度提高的同时，其质量和体积不断减小，已进入试验验证阶段，x射线脉冲星导航技术的发展前景十分可观。

导航精度是制约x射线脉冲星导航技术应用的瓶颈之一，导航精度取决于测量精度、模型精度、计时精度等方面。从验证导航方案可行性的角度出发，需要重点关注系统模型精度、测量信息精度对导航精度甚至稳定性的影响，集中体现在脉冲轮廓累积观测方法、轨道外推模型等方面。

调研表明，目前国内外脉冲星绝对导航算法方案中均采用了基于脉冲累积进行相位比较的方式产生导航测量。该方式首先在某个脉冲星观测时段内，通过航天器上的x射线探测器记录x射线脉冲光子到达时间(toa)，同时将通过轨道外推获得的位置信息用于建立时间转换方程，对光子toa观测量进行时间转换；其次，基于经过转换的光子toa数据，通过周期折叠方法获得观测脉冲轮廓；再次，比对观测脉冲轮廓和标准脉冲轮廓模版，得到的脉冲到达时间差(tdoa)；最后，结合航天器轨道动力学模型，通过滤波算法处理脉冲tdoa，获得航天器位置和速度的估计值。

上述基于脉冲累积进行相位比较的方式具有误差反馈直接的特点，针对长期项偏差具有较好的校正能力。然而，结合目前实际应用的情况，主要的问题在于：测量间隔长，脉冲累积固定，脉冲累积完成后才用于测量更新。这些因素造成的负面影响包括：一方面，目前的测量校正方式更新频率低，对于测量信息中随机误差的抑制能力较弱。因此，最终影响到定位精度的误差中，随机残留误差仍将保持在原始导航观测测量随机误差的水平。另一方面，由于测量间隔长，轨道外推误差将在测量间隔期间不断增大，导致光子toa转换出现偏差，进而引起脉冲轮廓累积误差，直接影响最终导航测量精度。

技术实现要素：

本发明的技术解决问题：克服现有技术的不足，提供一种基于动态脉冲累积窗口的x射线脉冲星导航方法，提高导航观测信息利用效率和测量更新频率，改善x射线脉冲星导航的精度和稳定性。

本发明的技术解决方案是：一种基于动态脉冲累积窗口的x射线脉冲星导航方法，包括：初始化阶段、状态预测阶段、光子测量处理阶段、脉冲轮廓动态累积阶段、tdoa求解阶段、滤波估计阶段。

所述初始化阶段步骤如下:

(1)设置光子测量参数，包括：当前累积初始时刻t0；发出光子的脉冲星的信号周期tp，简称为脉冲信号周期，tp为正数；每个脉冲周期内的时间窗格数nb，nb为正整数。设置轮廓累积参数，包括：累积脉冲周期数np，np为正整数。

(2)定义状态变量x：

x＝[rv]^t

其中r＝[rxryrz]和v＝[vxvyvz]分别为航天器在惯性坐标系中的位置矢量、速度矢量，rx表示r在惯性坐标系x轴上的分量，依此类推。

设置状态变量x的估计初值若已知航天器在惯性坐标系中的位置和速度先验值，则由航天器在惯性坐标系中位置先验估计值r0、速度先验估计值v0组成，否则可以取值为6维零向量；然后据此设置其他导航滤波参数，包括：对应的协方差矩阵p0，p0为6阶正定对角矩阵，每个对角元素可取值为中相应元素的估计方差；对应的系统噪声协方差矩阵q0，q0为6阶正定对角矩阵；对应的测量噪声协方差矩阵初值r0，r0为大于零的标量。

(3)设定当前滤波计算步数为k，k作为下标表示该变量在第k步的对应值，k的初值为1。记当前计算时刻为t。记滤波计算步长为tf，tf为正数。

(4)在t∈[t0,t0+np·tp]时间区间内，按顺序标记x射线脉冲星的脉冲信号周期，脉冲信号周期序号记为j，j＝1,2,...,np；同时，对每个脉冲信号周期，按顺序划分紧邻不重叠的时间窗格，每个时间窗格的宽度为tp/nb，时间窗格序号为i，i＝1,2,...,nb；定义nb×np维计数器数组c，c中的元素c(i,j)表示序号为j的脉冲周期中，序号为i的时间窗格内的光子累积计数，c(i,j)取值范围为非负整数，初值为0。

所述状态预测阶段步骤如下:

(5)根据第k-1步的状态估计值和误差方差阵pk-1，以及离散化卫星轨道动力学模型fk，预测当前，即第k步的状态变量和相应的误差方差阵pk|k-1，具体公式为：

其中，下标“k|k-1”表示根据第k-1步信息预测得到的第k步信息。若k＝1，和pk-1的初值分别为状态估计初值状态估计协方差矩阵初值p0。

离散化卫星轨道动力学模型fk按下式计算：

fk＝i6×6+dk·tf

离散状态转移矩阵dk为：

矩阵d21中的各元素分别为：

d51＝d42

d61＝d43

d62＝d53

其中μ为地球引力常数；j2为地球引力二阶带谐项系数，由外部给出；re为地球平均赤道半径，由外部给出；矩阵d21中的位置矢量r等于中的位置矢量部分。

所述光子测量处理阶段步骤如下：

对每一个到达航天器的光子，都执行步骤(a1)～步骤(a4)。

(a1)利用x射线导航敏感器toa测量值，记录第p个光子到达航天器的时间tp,sc，p为光子的序号，初值为0。

(a2)根据步骤(5)得出的和步骤(a1)获得的tp,sc，进行光子toa时间转换，得到的第p个光子到达太阳系质心的时间时间转换方程如下所示：

其中，的取值为即中的位置矢量部分；是矢量的大小；c表示真空光速；矢量n是发出光子的脉冲星视线方向单位矢量，由外部提供；d0是太阳系质心到脉冲星的距离，由外部提供；μs表示太阳引力场系数，由外部提供；b是太阳系质心相对于太阳系中心的位置矢量，由外部提供；b是b矢量的大小。

(a3)计算对应的归一化相位值mod(·)表示取余运算：

计算归一化相位对应的时间窗格序号ip，int(·)表示四舍五入取整运算：

以及当前光子对应的脉冲周期序号jp：

然后，设置计数器c(ip,jp)加1；

(a4)存储c(i,j),i＝1,2,...,nb,j＝1,2,...,np，k加1后继续执行步骤(a1)～(a4)，即对下一个观测到的光子进行处理。

所述脉冲轮廓动态累积阶段步骤如下：

(6)定义x射线脉冲星的观测累积时长toc为从t0开始到当前时刻t的时长，即

toc＝t-t0

若toc＝np·tp，则进入步骤(7)，否则返回步骤(5)。

(7)根据步骤(a4)的[t0,t0+np·tp]时间区间内的光子计数数据c(i,j)，得到观测脉冲速率函数由nb个按序排列的非负整数组成，计算公式如下：

所述tdoa求解阶段步骤如下：

(8)采用信号相位比较方法，计算观测脉冲速率函数与标准脉冲速率函数λ0(t)的相位差δφ。其中，标准脉冲速率函数表示t0为起始时刻的脉冲周期内太阳系质心处观测到的脉冲星标准计数序列，具体取值取决于所观测的脉冲星，通过外部提供，为已知量。要求标准脉冲速率函数的数据格式与观测脉冲速率函数的数据步长、数据长度保持一致，即由nb个按序排列的非负整数组成。

(9)计算观测脉冲速率函数与标准脉冲速率函数λ0(t)的对应的到达时间差tdoa：

tdoa＝δφ·tp

所述滤波估计阶段步骤如下：

(10)计算当前状态变量的估计值及对应的状态估计协方差矩阵pk，具体公式为：

其中，观测量yk为步骤(9)给出的tdoa值；kk表示滤波增益矩阵。测量函数的表达式为：

观测方程雅克比矩阵hk的计算方法为：

hk＝[h1101×3]

其中，

(11)将步骤(10)得到的状态变量估计值作为导航输出，同时判定若满足设定的计算停止条件，则结束计算，即完成导航计算。否则，设置窗口移动步长tm，tm为正数且大于滤波计算步长tf，用t0+tm替换t0，返回步骤(6)，根据窗口长度需要，重新设置累积周期数np。

本发明所提出的基于动态脉冲累积窗口的x射线脉冲星导航方法适用于基于脉冲累积进行相位比较的导航方式。与现有基于固定窗口脉冲累积的导航技术相比具有以下优点：

(1)针对现有技术中脉冲累积窗口设置不灵活的问题，本发明采用动态脉冲累积窗口，随着脉冲星观测的更新不断调整窗口位置，在每个窗口中都进行脉冲轮廓累积，可适应不同工况的需求；

(2)针对现有技术中无法有效抑制随机误差的问题，本发明通过累积窗口与滤波递推的同步移动，能够在每次滤波递推都获得测量信息，从而大大提高测量更新频率，提高导航测量带宽，有效抑制测量随机误差；

(3)针对现有技术中轨道外推误差影响过大的问题，本发明缩短了测量更新间隔，显著降低了测量更新间隔期间轨道外推误差的增大，等效于提高了导航系统模型的精度；

(4)在算法复杂度和计算量方面，本发明改进了现有技术中的测量数据利用方法，采用序贯方式处理连续的观测数据，在提高导航信息利用效率的同时，不增加单次量测更新所需的计算量和数据存储需求，有利于算法的快速自主实施；

(5)本发明继承现有技术中的光子观测和计数方法，不增加硬件设备和观测负担，具有很好的算法移植性和适应能力。

附图说明

图1为动态脉冲累积窗口时间轴示意图；

图2为基于动态脉冲累积窗口的导航流程图；

图3为导航测量tdoa输出时间历程图；

图4为导航误差时间历程图。

具体实施方式

本发明的基本思路为：提出一种基于动态脉冲累积窗口的x射线脉冲星导航方法，属于航天器制导、导航与控制领域，有别于目前x射线脉冲星导航方法中对x射线脉冲星光子信号进行固定窗口累积的方式，本发明采用动态设置的脉冲累积窗口，随着x射线光子脉冲采样的更新，可根据信号累积信噪比和观测更新频率的需要，与滤波递推计算的过程同步地调整累积窗口位置，同时更新脉冲轮廓累积数据，提高导航观测更新频率。本发明可以有效解决传统x射线脉冲星光子信号固定窗口累积方式中无法有效抑制随机误差、轨道外推误差影响过大、脉冲累积窗口设置不灵活等问题，提高导航精度和稳定性。

下面结合附图和具体实施例对本发明进行详细说明。

本发明提供一种基于动态脉冲累积窗口的x射线脉冲星导航方法，提高导航观测信息利用效率和测量更新频率，改善x射线脉冲星导航的精度和稳定性。所述方法依托的设备包括x射线导航敏感器、星载计算机以及其它辅助设备。主要计算流程如附图1，包括初始化阶段、光子测量处理阶段、脉冲轮廓动态累积阶段、tdoa求解阶段、滤波估计阶段等。具体实施方式如下：

(一)初始化阶段

(1)设置光子测量参数，包括：当前累积初始时刻t0；发出光子的脉冲星的信号周期tp，简称为脉冲信号周期，tp为正数；每个脉冲周期内的时间窗格数nb，nb为正整数。设置轮廓累积参数，包括：累积脉冲周期数np，np为正整数。

(2)定义状态变量x：

x＝[rv]^t

其中r＝[rxryrz]和v＝[vxvyvz]分别为航天器在惯性坐标系中的位置矢量、速度矢量，rx表示r在惯性坐标系x轴上的分量，依此类推。

设置状态变量x的估计初值若已知航天器在惯性坐标系中的位置和速度先验值，则由航天器在惯性坐标系中位置先验估计值r0、速度先验估计值v0组成，否则可以取值为6维零向量；然后据此设置其他导航滤波参数，包括：对应的协方差矩阵p0，p0为6阶正定对角矩阵，每个对角元素可取值为中相应元素的估计方差；对应的系统噪声协方差矩阵q0，q0为6阶正定对角矩阵；对应的测量噪声协方差矩阵初值r0，r0为大于零的标量。

(3)设定当前滤波计算步数为k，k作为下标表示变量在第k步的对应值，k的初值为1。记当前计算时刻为t。记滤波计算步长为tf，tf为正数。

(4)在t∈[t0,t0+np·tp]时间区间内，按顺序标记x射线脉冲星的脉冲信号周期，脉冲信号周期序号记为j，j＝1,2,...,np；同时，对每个脉冲信号周期，按顺序划分紧邻不重叠的时间窗格，每个时间窗格的宽度为tp/nb，时间窗格序号为i，i＝1,2,...,nb；定义nb×np维计数器数组c，c中的元素c(i,j)表示序号为j的脉冲周期中，序号为i的时间窗格内的光子累积计数，c(i,j)取值范围为非负整数，初值为0。

(二)脉冲轮廓动态累积阶段

(5)根据第k-1步的状态估计值和误差方差阵pk-1，以及离散化卫星轨道动力学模型fk，预测当前，即第k步的状态变量和相应的误差方差阵pk|k-1，具体公式为：

其中，下标“k|k-1”表示根据第k-1步信息预测得到的第k步信息。若k＝1，和pk-1的初值分别为状态估计初值状态估计协方差矩阵初值p0。

离散化卫星轨道动力学模型fk按下式计算：

fk＝i6×6+dk·tf

离散状态转移矩阵dk为：

矩阵d21中的各元素分别为：

d51＝d42

d61＝d43

d62＝d53

其中μ为地球引力常数；j2为地球引力二阶带谐项系数，由外部给出；re为地球平均赤道半径，由外部给出；矩阵d21中的位置矢量r等于中的位置矢量部分。

(三)光子测量处理阶段

对每一个到达航天器的光子，都执行步骤(a1)～步骤(a3)。

(a1)利用x射线导航敏感器toa测量值，记录第p个光子到达航天器的时间tp,sc，p为光子的序号，初值为0。

(a2)根据步骤(5)得出的和步骤(a1)获得的tp,sc，进行光子toa时间转换，得到的第p个光子到达太阳系质心的时间时间转换方程如下所示：

其中，的取值为即中的位置矢量部分；是矢量的大小；c表示真空光速；矢量n是发出光子的脉冲星视线方向单位矢量，由外部提供；d0是太阳系质心到脉冲星的距离，由外部提供；μs表示太阳引力场系数，由外部提供；b是太阳系质心相对于太阳系中心的位置矢量，由外部提供；b是b矢量的大小。

(a3)计算对应的归一化相位值mod(·)表示取余运算：

计算归一化相位对应的时间窗格序号ip，int(·)表示四舍五入取整运算：

以及当前光子对应的脉冲周期序号jp：

然后，设置计数器c(ip,jp)加1；

(a4)存储c(i,j),i＝1,2,...,nb,j＝1,2,...,np。k加1，继续执行步骤(a1)～(a4)，即对下一个观测到的光子进行处理。

(四)脉冲轮廓动态累积阶段

(6)定义x射线脉冲星的观测累积时长toc为从t0开始到当前时刻t的时长，即

toc＝t-t0

若toc＝np·tp，则进入步骤(7)，否则返回步骤(5)。

(7)根据步骤(a4)的[t0,t0+np·tp]时间区间内的光子计数数据c(i,j)，分别累加每个周期相同位置的时间窗格内的光子数，进行周期归一化，得到观测脉冲速率函数公式如下：

观测脉冲速率函数的物理含义是，以t0为起始时刻的脉冲周期内，对应每个ti时刻的光子计数c(ti,j)组成的序列，由nb个按序排列的非负整数组成。

(五)tdoa求解阶段

(8)采用信号相位比较方法，计算观测脉冲速率函数与标准脉冲速率函数λ0(t)的相位差δφ。其中，标准脉冲速率函数表示t0为起始时刻的脉冲周期内太阳系质心处观测到的脉冲星标准计数序列，具体取值取决于所观测的脉冲星，通过外部提供，为已知量。要求标准脉冲速率函数的数据格式与观测脉冲速率函数的数据步长、数据长度保持一致，即由nb个按序排列的非负整数组成。

(9)计算观测脉冲速率函数与标准脉冲速率函数λ0(t)的对应的到达时间差tdoa：

tdoa＝δφ·tp

(六)滤波估计阶段

(10)采用kalman滤波方法的算法框架，计算当前状态变量的估计值及对应的状态估计协方差矩阵pk，具体公式为：

其中，观测量yk为步骤(9)给出的tdoa值；kk表示滤波增益矩阵。测量函数的表达式为：

观测方程雅克比矩阵hk的计算方法为：

hk＝[h1101×3]

其中，

(11)将步骤(10)得到的状态变量估计值作为导航输出，同时判定若满足设定的计算停止条件，则结束计算，即完成导航计算。否则，设置窗口移动步长tm，tm为正数且大于滤波计算步长tf，用t0+tm替换t0，返回步骤(6)，根据窗口长度需要，重新设置累积周期数np。

如附图2，通过在步骤(11)对tm的按需设置，本发明上述算法可以通过累积窗口的动态移动，在选定的滤波计算步长中行脉冲轮廓累积，从而产生用于步骤(10)滤波估计过程的tdoa测量。由于窗口移动步长最短可等于测量更新步长，上述方式最大限度地提高了量测更新频率，有效抑制随机误差。同时，通过对np的设置，累积窗口的长度也可以根据需求来定义。足够长的累积窗口有利于保证累积信号的信噪比，提高tdoa测量精度，从而保证导航性能。

实施例

航天器轨道参照geo卫星数据，脉冲星选取、脉冲星导航敏感器参数等各特征参数如表1：

表1x射线脉冲星导航特征参数

按照表1的参数设置进行仿真，导航测量tdoa输出时间历程如附图3，其中选取仿真开始阶段1000s为例，对比了传统固定累积窗口方法和本发明动态累积窗口方法的导航测量输出标志，该标志为“1”表示有导航测量输出，为“0”表示没有导航测量输出。附图3上半部分为固定累积窗口方法，下半部分为动态累积窗口方法。可以明显看出，在经历前300秒左右的初始化阶段后，动态累积窗口方法的测量更新间隔即为窗口移动步长(10s)，与固定累积窗口方法约300s的测量更新周期相比，测量更新频率提高了约30倍。

导航误差时间历程如附图4，图中横坐标为仿真时间，历时5天。纵坐标为惯性系三轴方向的位置误差。可以看到，滤波估计约1个轨道周期后收敛，导航输出保持平稳，三轴位置误差大小不超过1km。图4表明，本发明通过较高的测量更新频率，使得导航输出保持平稳。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。