CN104296755B - 一种x射线脉冲星导航脉冲toa的确定方法 - Google Patents

Abstract

本发明涉及一种X射线脉冲星导航脉冲TOA的确定方法，该方法首先根据观测时间内得到的光子序列，在选取的折叠周期内进行脉冲轮廓折叠，并且进行归一化得到归一化强度矢量，再对该归一化强度矢量的协方差矩阵进行特征子空间分解得到噪声子空间矩阵，并根据标准脉冲轮廓构造搜索向量，然后由该噪声子空间矩阵和搜索矢量得到脉冲TOA的谱，并通过对脉冲TOA的谱进行谱峰值搜索得到脉冲TOA的值；该方法基于特征子空间分解进行脉冲TOA确定可以有效提高脉冲TOA的计算精度，并采用脉冲轮廓折叠提高脉冲轮廓估计值与标准脉冲轮廓的逼近程度，从而提高观测时间内的脉冲TOA计算精度。

Description

一种X射线脉冲星导航脉冲TOA的确定方法

技术领域

本发明涉及航天器自主导航技术领域，特别涉及一种X射线脉冲星导航脉冲TOA的确定方法。

背景技术

X射线脉冲星导航(XPNAV)是一种新概念的航天器自主导航技术，能够为近地轨道、深空和星际空间飞行的航天器提供位置、速度、时间和姿态等丰富的导航信息，具有重要的工程应用价值和战略研究意义，备受国际航天机构关注。从公开报道来看，X射线脉冲星导航最知名的计划是2004年美国国防先进技术发展局启动的“X射线导航与自主定位”(XNAV)研究计划，该计划目前已完成可行论证、关键技术攻关与地面验证，即将在国际空间站和高轨道卫星上开展空间飞行试验。其它公开的脉冲星导航计划还包括，2004年欧洲空间局(ESA)在ARIADNA空间技术预研计划支持下启动的“ESA深空探测器脉冲星导航研究计划”，进行了可行性研究；日本、俄罗斯和澳大利亚等国家也开展了相关的研究。

XPNAV的基本原理是将X射线脉冲星到达航天器的到达时间(TOA)作为基本观测量，利用建立在基准点(太阳系质心)的时间模型，计算同一脉冲到达基准点的TOA，利用脉冲TOA观测量与预测模型之间的差值，通过一定的导航算法，获得观测时刻航天器相对基准点的位置。因此只有进行精确的脉冲TOA确定，才能进行后续的航天器导航定位。

脉冲TOA是XPNAV的基本观测量，其确定精度是决定XPNAV定位、测速和定时精度的主要因素，如何高精度的求解脉冲TOA是XPNAV的一个十分重要的问题。Emadzadeh A A等人提出了一系列的脉冲TOA确定方法，在文献“Emadzadeh A A,Speyer J L.On modeling andpulse phase estimation of X-ray pulsars.IEEE Transactions on SignalProcessing,2010,58(9):4484–4495”和文献“Emadzadeh A A,Speyer J L.X-Ray Pulsar-Based Relative Navigation using Epoch Folding.IEEE Transactions on Aerospaceand Electronic Systems,2011:0018-9251”中提出将到达光子进行历元折叠得到折叠轮廓，前者利用非线性最小二乘估计而后者利用互相关技术实现折叠轮廓与标准轮廓的对比，得到脉冲TOA的值；文献“Emadzadeh A A,Speyer J L.Asymptotically efficientestimation of pulse time delay for X-ray pulsar based relative navigation.In:AIAAGN&C Conference,Chicago,IL,2009:1–12”中基于最大似然原理，提出直接利用光子数据通过最大化似然函数求得脉冲TOA的值。文献“Golshan A R,Sheikh S I.On pulsephase estimation and tracking of variable celestial X-ray sources.In:ION 63rdAnnual Meeting,Cambridge,MA,2007:413–422”提出利用栅格化搜索的数值算法来完成似然函数的最大化。Rinauro S等人在文献“Rinauro S,Colonnese S,ScaranoG.Fast near-maximum likelihood phase estimation of X-ray pulsars.Signal Processing,2013,93(1):326-331”中将脉冲TOA确定问题重构为一个循环移位参数估计问题，然后利用离散傅里叶变换完成似然函数的最大化，得到脉冲TOA的值。在国内，学者们也提出了一些高精度的脉冲TOA确定方法。谢振华等人在文献“谢振华，许录平，倪广仁.基于双谱的脉冲星累积脉冲轮廓时间延迟测量.物理学报,2008,57(10):6683-6688”中利用双谱抑制高斯噪声，在低信噪比下也具有较高的估计精度，但是该方法是针对射电脉冲星累积轮廓的脉冲TOA确定提出的，并不适用与XPNAV；苏哲等人在文献“苏哲，许录平，王婷等.一种新的脉冲星累积脉冲轮廓时间延迟测量算法.宇航学报,2011,32(6):1256-1261.”中提出了一种粗略估计和精确测量相结合的XPNAV脉冲TOA确定方法，通过抛物面内插法来提高精度。

但是，为了获得更好的导航效果，现有的脉冲TOA确定方法的精度仍有待进一步提高。

发明内容

本发明的目的在于克服现有技术的不足，提供一种X射线脉冲星导航脉冲TOA的确定方法，该方法基于特征子空间分解进行脉冲TOA的确定，可以有效提高脉冲TOA的计算精度，并采用脉冲轮廓折叠提高脉冲轮廓估计值与标准脉冲轮廓的逼近程度，从而提高观测时间内的脉冲TOA的计算精度。

本发明的上述目的是通过如下技术方案予以实现的：

一种X射线脉冲星导航脉冲TOA的确定方法，包括下列步骤：

(1)、根据设定的观测脉冲星序号，得到所述脉冲星的脉冲周期P和标准脉冲轮廓S(t)；

(2)、从观测时间起点τ₀开始，利用探测器接收光子，并记录所述光子的到达时间，在设定的观测时间τ_obs内，得到光子序列T_g＝[τ₁,τ₂,…,τ_k,…,τ_Ng]^T，其中τ_k为第k个光子到达探测器的时间，k＝1、2、…、Ng，Ng为观测过程中所述探测器接收到的光子总数；

(3)、根据步骤(2)得到的光子序列，进行光子历元折叠并进行归一化，得到折叠轮廓归一化强度矢量R＝[r₁,r₂,…,r_N]^T，具体处理过程如下：

(3a)、在观测时间段τ₀～τ₀+τ_obs内选取折叠起点为τ_0f，并以脉冲周期P为折叠周期，得到折叠时间段τ_0f～τ_0f+P，其中τ₀≤τ_0f≤τ₀+τ_obs-P；

(3b)、将步骤(3a)得到的折叠时间段均匀划分为N个折叠时间格，其中，第n个折叠时间格为τ_0f+(n-1)τ_b～τ_0f+nτ_b，其中，n＝1，2，…,N，所述折叠时间格长度

(3c)、在k＝1、2、…、Ng时，根据如下判断方法确定各光子到达时间在折叠时间段内的位置，并统计得到各折叠时间格内的光子数：

如果第k个光子到达探测器的时间τ_k满足条件：则判断所述第k个光子到达时间τ_k位于第n个折叠时间格内，并将所述第n个折叠时间格内的光子数C(n)加1；其中表示计算结果的小数部分的数值，C(n)的初始值为0；

(3d)、对步骤(3c)中统计得到的各折叠时间格内的光子数进行归一化，得到折叠轮廓归一化强度矢量R＝[r₁,r₂,…,r_N]^T，其中：

其中，C_min为C(1)、C(2)、…、C(N)内的最小值；

(4)、根据步骤(3)得到的折叠轮廓归一化强度矢量R＝[r₁,r₂,…,r_N]^T，得到协方差矩阵并对所述协方差矩阵进行特征值分解得到N个特征值，以及所述N个特征值对应的N个特征矢量；其中，所述N个特征值由大到小依次为λ₁、λ₂、…、λ_N，即λ₁>λ₂>…>λ_N；与特征值λ₁、λ₂、…、λ_N对应的特征向量依次为u₁、u₂、…、u_N；

(5)、根据步骤(4)中得到的特征向量u₂、u₃、…、u_N，得到噪声子空间矩阵U＝[u₂,u₃,…,u_N]；

(6)、根据步骤(1)得到的标准脉冲轮廓S(t)构造搜索向量A(φ₁)、A(φ₂)、…、A(φ_n)、…、A(φ_N)，其中：

A(φ_n)＝[s_n,s_n+1,…,s_N,s₁,s₂,…,s_n-1]^T

其中s₁～s_N为对标准脉冲轮廓S(t)进行N点均匀采样得到的N点标准脉冲轮廓值；n＝1，2，…,N；

(7)、根据步骤(5)计算得到的噪声子空间矩阵U＝[u₂,u₃,…,u_N]和步骤(6)构造得到的搜索向量A(φ₁)、A(φ₂)、…、A(φ_n)、…、A(φ_N)，计算脉冲TOA在不同搜索向量下的谱f(φ₁)、f(φ₂)、…、f(φ_n)、…、f(φ_N)，其中

(8)、在f(φ₁)～f(φ_N)中查找最大值f_max，如果f_max＝f(φ_j)，则导航脉冲TOA为P·φ_j，其中j＝1,2,…,或N。

在上述的X射线脉冲星导航脉冲TOA的确定方法中，在步骤(3c)中经过光子历折叠统计得到的每个折叠时间格内的光子数C(1)、C(2)、…、C(N)为折叠得到脉冲轮廓。

在上述的X射线脉冲星导航脉冲TOA的确定方法中，观测时间τ_obs与X射线脉冲星导航脉冲TOA的计算精度成正比，即观测时间越长，所述导航脉冲TOA的计算精度越高。

本发明与现有技术相比具有如下有益效果：

(1)、本发明首先计算得到折叠脉冲轮廓归一化强度矢量，并对该矢量的协方差矩阵进行特征子空间分解得到噪声矢量，由此得到脉冲TOA的谱，通过谱搜索得到脉冲TOA的值，可以显著提高X射线脉冲星的脉冲TOA确定精度；

(2)、本发明在观测时间内设定折叠周期时段，将各观测时间内的光子折叠到折叠周期内，由此得到脉冲轮廓的估计值，即折叠脉冲轮廓，可以有效提高观测时间内的脉冲轮廓确定精度，因而提高脉冲TOA的确定精度。

附图说明

图1为本发明基于特征子空间分解的X射线脉冲星脉冲TOA确定方法的流程图；

图2为本发明实施例得到的X射线脉冲星脉冲相位谱图；

图3为现有方法1得到的X射线脉冲星脉冲相位谱图；

图4为本发明与现有方法1在不同观测时间长度内进行脉冲TOA确定的精度对比仿真结果。

具体实施方式

下面结合附图和具体实施例对本发明作进一步的详细描述：

本发明提供的X射线脉冲星导航脉冲TOA的确定方法，可以满足X射线脉冲星导航的精度指标要求，从而实现自主导航。

由图1所示的本发明方法流程图可知，本发明的射线脉冲星导航脉冲TOA的确定方法包括下列步骤：

(1)、由于每个脉冲星都对应确定的脉冲周期和标准脉冲轮廓，因此可以根据设定的观测脉冲星序号得到所述脉冲星的脉冲周期P和标准脉冲轮廓S(t)，用于以下步骤中脉冲轮廓折叠和搜索矢量确定；

该实施例中，选取的脉冲星为PSRB1744-24A，其脉冲周期P＝0.0116s；

(2)、从观测时间起点τ₀开始，利用探测器接收光子，并记录所述光子的到达时间，在设定的观测时间τ_obs内，得到光子序列T_g＝[τ₁,τ₂,…,τ_k,…,τ_Ng]^T，其中τ_k为第k个光子到达探测器的时间，k＝1、2、…、Ng，Ng为观测过程中所述探测器接收到的光子总数；

在进行脉冲TOA计算过程中，观测时间τ_obs与X射线脉冲星导航脉冲TOA的计算精度成正比，即观测时间越长，所述导航脉冲TOA的计算精度越高，在该实施例中，设定的观测时间观测时间τ_obs为2s；

(3)、根据步骤(2)得到的光子序列，进行光子历元折叠并进行归一化，得到折叠轮廓归一化强度矢量R＝[r₁,r₂,…,r_N]^T，具体处理过程如下：

(3a)、在观测时间段τ₀～τ₀+τ_obs内选取折叠起点为τ_0f，并以脉冲周期P为折叠周期，得到折叠时间段τ_0f～τ_0f+P，其中τ₀≤τ_0f≤τ₀+τ_obs-P；

(3b)、将步骤(3a)得到的折叠时间段均匀划分为N个折叠时间格，其中，第n个折叠时间格为τ_0f+(n-1)τ_b～τ_0f+nτ_b，其中，n＝1，2，…,N，所述折叠时间格长度

在该实施例中，折叠时间格数N设定为1024；

(3c)、在k＝1、2、…、Ng时，根据如下判断方法确定各光子到达时间在折叠时间段内的位置，并统计得到各折叠时间格内的光子数：

如果第k个光子到达探测器的时间τ_k满足条件：则判断所述第k个光子到达时间τ_k位于第n个折叠时间格内，并将所述第n个折叠时间格内的光子数C(n)加1；其中表示计算结果的小数部分的数值，C(n)的初始值为0；

经过以上光子历折叠处理得到的每个折叠时间格内的光子数C(1)、C(2)、…、C(N)为折叠得到脉冲轮廓；

(3d)、对步骤(3c)中统计得到的各折叠时间格内的光子数进行如下的归一化处理，得到折叠轮廓归一化强度矢量R＝[r₁,r₂,…,r_N]^T，其中：

其中，C_min为C(1)、C(2)、…、C(N)内的最小值；

通过以上完成折叠轮廓归一化后，归一化脉冲相位为该相位乘以脉冲周期P可以得到脉冲TOA，因此进行脉冲TOA的确定与脉冲相位的确定是等效的。

(4)、根据步骤(3)得到的折叠轮廓归一化强度矢量R＝[r₁,r₂,…,r_N]^T，得到协方差矩阵并对所述协方差矩阵进行特征值分解得到N个特征值，以及所述N个特征值对应的N个特征矢量；其中，所述N个特征值由小到大依次为λ₁、λ₂、…、λ_N，即λ₁>λ₂>…>λ_N；与特征值λ₁、λ₂、…、λ_N对应的特征向量依次为u₁、u₂、…、u_N；

(5)、在步骤(4)的特征分解结果中，最大的特征值λ₁为脉冲星特征值，与λ₁对应的特征矢量u₁可以张成脉冲星子空间，而剩下的N-1个特征值λ₂、…、λ_N为噪声特征值，因此与其对应的特征矢量u₂、…、u_N可以张成噪声子空间；根据步骤(4)中得到的特征向量u₂、u₃、…、u_N，得到噪声子空间矩阵U＝[u₂,u₃,…,u_N]；

(6)、根据步骤(1)得到的标准脉冲轮廓S(t)构造搜索向量A(φ₁)、A(φ₂)、…、A(φ_n)、…、A(φ_N)，其中：

A(φ_n)＝[s_n,s_n+1,…,s_N,s₁,s₂,…,s_n-1]^T

其中s₁～s_N为对标准脉冲轮廓S(t)进行N点均匀采样得到的N点标准脉冲轮廓值；n＝1，2，…,N；

(7)、根据步骤(5)计算得到的噪声子空间矩阵U＝[u₂,u₃,…,u_N]和步骤(6)构造得到的搜索向量A(φ₁)、A(φ₂)、…、A(φ_n)、…、A(φ_N)，计算脉冲TOA的谱f(φ₁)、f(φ₂)、…、f(φ_n)、…、f(φ_N)，其中n＝1，2，…,N；

(8)、在f(φ₁)～f(φ_N)中查找最大值f_max，如果f_max＝f(φ_j)，则导航脉冲TOA为P·φ_j，其中归一化的脉冲相位值j＝1,2,…,或N。

根据以上所述的方法实现步骤进行脉冲TOA确定的仿真分析，得到本发明的脉冲相位的谱图如图2所示，其中，谱峰值对应的归一化脉冲相位φ＝0.3999，即脉冲TOA的值为P·φ＝0.00463884s。

然后利用文献“Emadzadeh A A,Speyer J L.X-Ray Pulsar-Based RelativeNavigation using Epoch Folding.IEEE Transactions on Aerospace and ElectronicSystems,2011:0018-9251”提出的现有方法1，即采用互相关法进行脉冲TOA的确定，对该方法进行仿真得到的脉冲相位的谱如图3所示，其中，该谱峰值对应的归一化脉冲相位为φ＝0.3950，即脉冲TOA的值为P·φ＝0.004582s。

通过以上对本发明方法和现有方法1的仿真结果进行对比可以，在观测时间长度为2s时，两种方法均能正确确定出脉冲TOA，但本发明方法的脉冲相位确定谱较为尖锐，因此结果更准确。

以下观测时间长度分别取0.1s,0.2s,0.5s,1s,2s,5s,10s,20s,50s或100s，并在以上观测时间长度取值条件下分别进行500次蒙特卡洛试验，得到本发明方法和现有方法的脉冲TOA确定误差的比较结果，由图4所示的仿真对比结果可以看出，本发明可以有效提高脉冲TOA的确定精度，并且在较短的观测时间下具有较好的结果。

以上所述，仅为本发明最佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

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

Claims (3)

1.一种X射线脉冲星导航脉冲TOA的确定方法，其特征在于包括下列步骤：

(1)、根据设定的观测脉冲星序号，得到所述脉冲星的脉冲周期P和标准脉冲轮廓S(t)；

(2)、从观测时间起点τ₀开始，利用探测器接收光子，并记录所述光子的到达时间，在设定的观测时间τ_obs内，得到光子序列T_g＝[τ₁,τ₂,…,τ_k,…,τ_Ng]^T，其中τ_k为第k个光子到达探测器的时间，k＝1、2、…、Ng，Ng为观测过程中所述探测器接收到的光子总数；

(3)、根据步骤(2)得到的光子序列，进行光子历元折叠并进行归一化，得到折叠轮廓的归一化强度矢量R＝[r₁,r₂,…,r_N]^T，具体处理过程如下：

(3a)、在观测时间段τ₀～τ₀+τ_obs内选取折叠起点为τ_0f，并以脉冲周期P为折叠周期，得到折叠时间段τ_0f～τ_0f+P，其中τ₀≤τ_0f≤τ₀+τ_obs-P；

(3b)、将步骤(3a)得到的折叠时间段均匀划分为N个折叠时间格，其中，第n个折叠时间格为τ_0f+(n-1)τ_b～τ_0f+nτ_b，其中，n＝1，2，…,N，所述折叠时间格长度

(3c)、在k＝1、2、…、Ng时，根据如下判断方法确定各光子到达时间在折叠时间段内的位置，并统计得到各折叠时间格内的光子数：

如果第k个光子到达探测器的时间τ_k满足条件：则判断所述第k个光子到达时间τ_k位于第n个折叠时间格内，并将所述第n个折叠时间格内的光子数C(n)加1；其中表示计算结果的小数部分的数值，C(n)的初始值为0；

(3d)、对步骤(3c)中统计得到的各折叠时间格内的光子数进行如下的归一化处理，得到折叠轮廓归一化强度矢量R＝[r₁,r₂,…,r_N]^T，其中：

$r_n=\frac{P}{{\mathrm{\τ}}_b}\frac{C\left(n\right)-C_{min}}{\underset{m=1}{\overset{N}{\Σ}}\[C\left(m\right)-C_{min}\]},n=1,2,\mathrm{...},N$

其中，C_min为C(1)、C(2)、…、C(N)内的最小值；

(4)、根据步骤(3)得到的折叠轮廓归一化强度矢量R＝[r₁,r₂,…,r_N]^T，计算协方差矩阵并对所述协方差矩阵进行特征值分解得到N个特征值，以及所述N个特征值对应的N个特征矢量；其中，所述N个特征值由大到小依次为λ₁、λ₂、…、λ_N，即λ₁>λ₂>…>λ_N；与特征值λ₁、λ₂、…、λ_N对应的特征向量依次为u₁、u₂、…、u_N；

(5)、根据步骤(4)中得到的特征向量u₂、u₃、…、u_N，得到噪声子空间矩阵U＝[u₂,u₃,…,u_N]；

(6)、根据步骤(1)得到的标准脉冲轮廓S(t)构造搜索向量A(φ₁)、A(φ₂)、…、A(φ_n)、…、A(φ_N)，其中：

A(φ_n)＝[s_n,s_n+1,…,s_N,s₁,s₂,…,s_n-1]^T

其中s₁～s_N为对标准脉冲轮廓S(t)进行N点均匀采样得到的N点标准脉冲轮廓值；n＝1，2，…,N；

(7)、根据步骤(5)计算得到的噪声子空间矩阵U＝[u₂,u₃,…,u_N]和步骤(6)构造得到的搜索向量A(φ₁)、A(φ₂)、…、A(φ_n)、…、A(φ_N)，计算脉冲TOA估计谱f(φ₁)、f(φ₂)、…、f(φ_n)、…、f(φ_N)，其中n＝1，2，…,N；

(8)、在f(φ₁)～f(φ_N)中查找最大值f_max，如果f_max＝f(φ_j)，则导航脉冲TOA为P·φ_j，其中j＝1,2,…,或N。

2.根据权利要求1所述的一种X射线脉冲星导航脉冲TOA的确定方法，其特征在于，在步骤(3c)中经过光子历元折叠统计得到的每个折叠时间格内的光子数C(1)、C(2)、…、C(N)为折叠得到脉冲轮廓。

3.根据权利要求1所述的一种X射线脉冲星导航脉冲TOA的确定方法，其特征在于，观测时间τ_obs与X射线脉冲星导航脉冲TOA的计算精度成正比，即观测时间越长，所述导航脉冲TOA的计算精度越高。