[发明专利]一种基于动态脉冲累积窗口的X射线脉冲星导航方法

权利要求书

1.一种基于动态脉冲累积窗口的X射线脉冲星导航方法，其特征在于包括：初始化阶段、状态预测阶段、光子测量处理阶段、脉冲轮廓动态累积阶段、TDOA求解阶段、滤波估计阶段，

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

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

(2)定义状态变量x：

x＝[r v]^T

其中r＝[r_x r_y r_z]和v＝[v_x v_y v_z]分别为航天器在惯性坐标系中的位置矢量、速度矢量，r_x表示r在惯性坐标系x轴上的分量，r_y表示r在惯性坐标系y轴的分量，r_z表示r在惯性坐标系z轴的分量；

设置状态变量x的估计初值包括航天器在惯性坐标系中位置先验估计值r₀、速度先验估计值v₀；然后设置其他导航滤波参数，包括：对应的协方差矩阵P₀，P₀为6阶正定对角矩阵；对应的系统噪声协方差矩阵Q₀，Q₀为6阶正定对角矩阵；对应的测量噪声协方差矩阵初值R₀，R₀为大于零的标量；

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

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

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

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

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

计算离散化卫星轨道动力学模型F_k；

离散状态转移矩阵D_k为：

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

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

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

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

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

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

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

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

然后，设置计数器C(i_p,j_p)加1；C(i_p,j_p)的初值为0

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

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

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

t_oc＝t-t₀

若t_oc＝N_p·T_p，则进入步骤(7)，否则返回步骤(5)；

(7)根据步骤(a4)的[t₀,t₀+N_p·T_p]时间区间内的光子计数数据C(i,j)，得到观测脉冲速率函数由N_b个按序排列的非负整数组成；

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

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

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

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

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

其中，观测量y_k为步骤(9)给出的TDOA值；K_k表示滤波增益矩阵；为测量函数；

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

H_k＝[H₁₁ 0_1×3]

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