CN104567937A - 一种高精度的x射线脉冲星信号模拟方法 - Google Patents

Abstract

本发明公开了一种高精度的X射线脉冲星信号模拟方法，该方法首先确定模拟的起始时刻和结束时刻，并将起始时刻视为一个光子达到时刻，根据脉冲星的相位预测模型计算得到光子相位，再利用反函数法递推得到下一个光子的相位，并将该相位代入脉冲星的相位预测模型中得到方程，求解该方程得到下一个光子的到达时刻，直到模拟过程结束；发明模拟产生的X射线脉冲星光子信号具有脉冲星自转频率缓慢变化的特性，能够用于构建高精度的脉冲星导航地面试验系统和脉冲星周期搜索、轮廓解算、TOA估计和长时间自主导航算法的验证。

Description

一种高精度的X射线脉冲星信号模拟方法

技术领域

本发明涉及一种高精度的X射线脉冲星信号模拟方法，属于信号模拟和信号处理领域。

背景技术

X射线脉冲星导航(XPNAV)是一种全新的自主导航方式，具有隐蔽性好、安全性高、抗干扰能力强，导航精度高等诸多优点，能够为近地轨道、深空和星际空间飞行的航天器提供位置、速度、时间和姿态等丰富的导航信息，具有重要的工程应用价值和战略研究意义，备受国际航天机构关注。2004年，美国国防先进技术发展局启动的“X射线导航与自主定位”(XNAV)研究计划，开展脉冲星仿真方法研究和原理样机研制等工作，该计划目前已完成可行论证、关键技术攻关与地面验证，即将在国际空间站和高轨道卫星上开展空间飞行试验。由于X射线无法穿透稠密的大气层，地面实验室无法观测到X射线脉冲星的辐射信号，因此卫星搭载探测器对X射线脉冲星进行观测。然而搭载试验具有复杂度高，耗时长，技术难度大等缺点，不适用于X射线脉冲星导航系统的研究和算法的验证。因此，目前的研究工作主要采用地面模拟产生的X射线脉冲星信号。

目前涉及X射线脉冲信号模拟方法的文献主要有：(1)“一种基于随机单光子发射机制的X射线脉冲星信号模拟源”(专利号：CN103697908A)，(2)“多脉冲星信号模拟器”(申请号：200910023707.3)，(3)“一种脉冲星频率信号模拟器”(专利号：CN102759884A)，(4)“刘利,郑伟,汤国建,孙守明.基于X射线脉冲星的导航半实物仿真系统[J].国防科技大学学报.2012,34(5)”，(5)“孙海峰,谢楷,李小平等.高稳定度X射线脉冲星信号模拟[J].物理学报.2013(10):518-528”，(6)“周峰,吴光敏,赵宝升等.基于X射线脉冲星导航的模拟调制仿真源研究[J].物理学报,2013(11)”，(7)“Emadzadeh A,Speyer J L.On modeling and pulse phase estimation of X-raypulsars.IEEE Transactions on Signal Processing,2010,58(9):4484–4495”，(8)“ZHANG Hua,XU Luping,SONG Shibin,Jiao Rong[J].Journal of Physics.98(2014)189–200.”，(9)“J.Sala,A.Urruela,X.Villares,R.Estalella,Feasibility Study for a Spacecraft Navigation System Relying On Pulsar TimingInformation,European Space Agency Advanced Concepts Team ARIADNA Study,2004.”以上文献设计的脉冲星信号模拟系统或模拟方法都未考虑脉冲星自转频率存在缓慢变化的特点。只有文献“苏哲.X射线脉冲星导航信号处理方法和仿真实验系统研究[D].2011.”考虑到脉冲星自转频率存在变化的特点，然而文献中还是假设在小时间段内脉冲星的频率是不变的，未提出针对频率缓慢变化特性的有效模拟方法。

脉冲星自主导航时间跨度大，频率的微小误差也会导致TOA巨大的偏差，因此只有准确地模拟脉冲星的频率变化特性，才能为脉冲星导航系统的研制和导航算法的验证提供真实有效的脉冲星信号。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供了一种高精度的X射线脉冲星信号模拟方法，模拟产生的X射线脉冲星光子信号具有脉冲星自转频率缓慢变化的特性，能够用于构建高精度的脉冲星导航地面试验系统和脉冲星周期搜索、轮廓解算、TOA估计和长时间自主导航算法的验证。

本发明的技术解决方案是：

一种高精度的X射线脉冲星信号模拟方法，步骤如下：

步骤S1：确定模拟的起始时刻t₀，模拟结束时刻t_f，记t_n为第n个脉冲星光子信号的到达太阳系质心SSB处的时间，令n的初始值为0；

步骤S2：根据脉冲星相位预测模型计算得到t_n时刻对应的脉冲星相位φ(t_n)；

步骤S3：由步骤S2得到的t_n时刻对应的脉冲星相位φ(t_n)，采用反函数法递推得到第n+1个脉冲星光子信号到达SSB处时对应的脉冲星相位φ(t_n+1)；

步骤S4：将步骤S3中得到的第n+1脉冲星光子信号对应的脉冲星相位φ(t_n+1)，代入脉冲星相位预测模型得到方程，求解该方程得到第n+1个脉冲星光子信号到达SSB处的时间t_n+1；

步骤S5：n的值加1，判断t_n是否小于模拟结束时刻t_f，当t_n<t_f时，转步骤S2；当t_n≥t_f时，结束本次模拟过程，此时得到脉冲星光子信号时间序列{t_0,t₁,t₂,…t_n,t_n+1,....t_f}，即实现了X射线脉冲星信号模拟。

脉冲星相位预测模型为

$\mathrm{\&phi;}\left(t_{n+1}\right)=\mathrm{\&phi;}\left(t_{\mathrm{ref}}\right)+f[t_{n+1}-t_{\mathrm{ref}}]+\underset{m=2}{\overset{m=N}{\mathrm{\&Sigma;}}}\frac{\frac{d^{m-1}f}{{\mathrm{dt}}^{m-1}}}{m!}{(t_{n+1}-t_{\mathrm{ref}})}^m$

其中，N是脉冲星相位预测模型的阶次，N取值范围为[2，+∞)，f和分别是脉冲星的自转频率和自转频率的m-1阶倒数值，φ(t)是t时刻脉冲到达SSB处的相位，φ(t_ref)是t_ref参考时刻脉冲到达SSB处的相位，t、t_ref和t_n都采用太阳系质心力学时TDB作为时间系统。

反函数递推公式如下：

φ(t_n+1)＝Λ^-1(Λ(φ(t_n))-lnU)

其中，U为[0,1]区间上均匀分布随机变量，Λ(i)为X射线脉冲星的累积泊松速率函数，Λ^-1(i)为累积泊松速率函数的反函数

本发明与现有技术相比的有益效果是：

(1)本发明得到的每一个脉冲星光子信号都和脉冲星的相位预测模型相关，整个脉冲星光子到达时间序列就包含了脉冲星自转频率缓慢变化的特性，符合脉冲星的实际情况。

(2)本发明模拟的X射线脉冲星光子达到时间序列的周期和脉冲星的实际自转周期误差很小，具有模拟精度高的优点。

(3)本发明光子通过递推的方法逐个产生，模拟的过程精确可控。

(4)本发明将反函数公式中的累积泊松速率函数Λ(i)替换成不同脉冲星的累积泊松速率函数，将相位预测模拟替换成相应的脉冲星相位预测模型，即可模拟不同脉冲星的光子信号，具有通用性强的优点。

附图说明

图1为本发明的算法步骤流程图；

图2为本发明反函数法光子相位递推算法的求解示意图；

具体实施方式

如图1所示是本发明方法实施步骤，本发明提供了一种高精度的X射线脉冲星信号模拟方法，步骤如下：

步骤S1：确定模拟的起始时刻t₀，模拟结束时刻t_f，记t_n为第n个脉冲星光子信号的到达太阳系质心SSB处的时间，令n的初始值为0；

步骤S2：根据脉冲星相位预测模型计算得到t_n时刻对应的脉冲星相位φ(t_n)；

脉冲星相位预测模型为

$\mathrm{\&phi;}\left(t_{n+1}\right)=\mathrm{\&phi;}\left(t_{\mathrm{ref}}\right)+f[t-t_{\mathrm{ref}}]+\underset{m=2}{\overset{m=N}{\mathrm{\&Sigma;}}}\frac{\frac{d^{m-1}f}{{\mathrm{dt}}^{m-1}}}{m!}{(t-t_{\mathrm{ref}})}^m$

其中，N是脉冲星相位预测模型的阶次，N取值范围为[2，+∞)，模型阶次N的选择随脉冲星自转频率的稳定度和计算精度的需要选择合适的值，对于大多数的应用N取3就能达到足够高的精度，f和分别是脉冲星的自转频率和自转频率的各阶倒数值，查找脉冲星星历表可以得到，φ(t)是t时刻脉冲到达SSB处的相位，φ(t_ref)是t_ref参考时刻脉冲到达SSB处的相位，t、t_ref和t_n都采用太阳系质心力学时TDB作为时间系统；

步骤S3：由步骤S2得到的t_n时刻对应的脉冲星相位φ(t_n)，采用反函数法递推得到第n+1个脉冲星光子信号到达SSB处时对应的脉冲星相位φ(t_n+1)；

反函数递推公式如下：

φ(t_n+1)＝Λ^-1(Λ(φ(t_n))-lnU)

其中，U为[0,1]区间上均匀分布随机变量，Λ(i)为X射线脉冲星的累积泊松速率函数，Λ^-1(i)为累积泊松速率函数的反函数；

由于脉冲星的泊松速率函数很难用函数解析表达式描述，同样地累积泊松速率函数和累积泊松速率函数的反函数也无法得到函数解析表达式，因此采用数值法求解。

步骤S4：将步骤S3中得到的第n+1脉冲星光子信号对应的脉冲星相位φ(t_n+1)，代入脉冲星相位预测模型得到方程：

$\mathrm{\&phi;}\left(t_{n+1}\right)=\mathrm{\&phi;}\left(t_{\mathrm{ref}}\right)+f[t_{n+1}-t_{\mathrm{ref}}]+\underset{m=2}{\overset{m=N}{\mathrm{\&Sigma;}}}\frac{\frac{d^{m-1}f}{{\mathrm{dt}}^{m-1}}}{m!}{(t_{n+1}-t_{\mathrm{ref}})}^m$

求解该方程，得到第n+1个脉冲星光子信号到达SSB处的时间t_n+1；

当方程的阶次N等于2时，此时方程是一个二次方程，可以用二次方程求根公式求解；当方程阶数N大于2时，可以用迭代法求解，由于相位演化具有单调特性，因此相位预测模型函数是一个光滑的单调函数，以任意时刻作为迭代起始点总能收敛到真值。

步骤S5：n的值加1，判断t_n是否小于模拟结束时刻t_f，当t_n<t_f时，转步骤S2；当t_n≧t_f时，结束本次模拟过程，此时得到脉冲星光子信号时间序列{t_0,t₁,t₂,…t_n,t_n+1,....t_f}；从而实现了高精度的X射线脉冲星信号模拟。

下面结合附图并举实施例，对本发明进行具体描述。

步骤S1：本实施例以脉冲星PSR B0531+21为例，模拟起始时间t₀为49128.000000191(MJD)，模拟时长180s，对应模拟结束时刻t_f为49128.002083524(MJD)，记t_n为第n个脉冲星光子信号的到达SSB(太阳系质心)处的时间，令n的初始值为0；

步骤S2：根据第n个光子的到达时刻t_n，查找脉冲星星历表，得到适用于当前时刻的脉冲星相位预测模型参数，对于脉冲星PSR B0531+21，N取3时就能达到足够高的精度，此时得到的脉冲星PSR B053+21相位预测模型为：

$\mathrm{\&phi;}\left(t\right)=\mathrm{\&phi;}\left(t_{\mathrm{ref}}\right)+f[t-t_{\mathrm{ref}}]+\frac{\stackrel{\&CenterDot;}{f}}{2}{[t-t_{\mathrm{ref}}]}^2+\frac{\stackrel{\&CenterDot;\&CenterDot;}{f}}{6}{[t-t_{\mathrm{ref}}]}^3---\left(1\right)$

将t_n代入(1)式，计算得到t_n时刻对应的脉冲星相位φ(t_n)，其中f是脉冲星的自转频率，是脉冲星自转频率一阶项，脉冲星自转频率二阶项，φ(t_ref)是t_ref参考时刻脉冲到达SSB处的相位，f、和t_ref通过查找脉冲星星历表可以得到，t_ref和t_n都采用太阳系质心力学时TDB作为时间系统；

步骤S3：由步骤S2得到的t_n时刻对应的脉冲星相位φ(t_n)，采用反函数法递推得到第n+1个脉冲星光子信号到达SSB处时对应的脉冲星相位φ(t_n+1)，递推公式如下：

φ(t_n+1)＝Λ^-1(Λ(φ(t_n))-lnU) (2)

其中，U为[0,1]区间上均匀分布随机变量，Λ(i)为X射线脉冲星的累积泊松速率函数，Λ^-1(i)为累积泊松速率函数的反函数。

由于脉冲星的泊松速率函数很难用函数解析表达式描述，同样地累积泊松速率函数和累积泊松速率函数的反函数也无法得到函数解析表达式，因此采用数值法求解。图2是反函数法光子相位递推算法的求解示意图，首先根据脉冲星的泊松速率得到相位—泊松速率值对应关系即Λ(i)，采用求和运算代替积分运算得到相位—累积泊松速率值对应关系即Λ^-1(i)，再根据前一个光子的相位φ(t_n)，查找相位—泊松速率值对应关系，得到Λ(φ(t_n))的结果，接着产生(0,1)上的均匀随机数U，计算得到下一个光子的对应的累积泊松速率值Λ(φ(t_n))-lnU，最后查找相位—累积泊松速率值对应关系，求解得到Λ(i)递推得到Λ^-1(Λ(φ(t_n))-lnU)的结果即下一个光子的相位φ(t_n+1)。

步骤S4：将步骤S3中得到的第n+1脉冲星光子信号对应的脉冲星相位φ(t_n+1)，代入(1)式的左边，得到方程

$\mathrm{\&phi;}\left(t_{n+1}\right)=\mathrm{\&phi;}\left(t_0\right)+f[t-t_0]+\frac{\stackrel{\&CenterDot;}{f}}{2}{{[t}_{n+1}-t_0]}^2+\frac{\stackrel{\&CenterDot;\&CenterDot;}{f}}{6}{[t_{n+1}-t_0]}^3---\left(3\right)$

求解方程(3)，得到第n+1个脉冲星光子信号到达SSB处的时间t_n+1；

步骤S5：n的值加1，判断t_n是否小于模拟结束时刻t_f，当t_n<t_f时，转步骤S2；当t_n≧t_f时，结束本次模拟过程，此时得到脉冲星光子信号时间序列{t_0,t₁,t₂,…t_n,t_n+1,....t_f}；

采用以上步骤，得到脉冲星PSR B0531+21起始时刻为50200.000000251、50645.000000104、50849.000000344、51224.000000271，持续时间为180s的4组脉冲星模拟光子到达时间序列。采用周期搜索算法，搜素模拟光子时间序列的周期，结果如表1所示：

表1.模拟光子信号时间序列的周期

仿真起始时间(MJD)	星历表的周期(f0)	模拟光子序列的周期
			49128.000000191	29.9245763662884	29.92458
50200.000000251	29.8897220169126	29.88972
			50645.000000104	29.8752812172442	29.87528
50849.000000344	29.8686661308798	29.86867
			51224.000000271	29.8565144334781	29.85651

由表可知模拟光子序列的平均周期和脉冲星星历表的周期很接近，由于脉冲星周期是持续缓慢变化的，而对模拟光子序列进行周期搜索得到的周期是整个时间序列的平均周期，因此和脉冲星星历表中的周期存在偏差。从结果可以看出模拟的光子序列具有脉冲星周期缓慢变化的特性，符合脉冲星的真实情况。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

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

Claims (3)

1.一种高精度的X射线脉冲星信号模拟方法，其特征在于步骤如下：

步骤S1：确定模拟的起始时刻t₀，模拟结束时刻t_f，记t_n为第n个脉冲星光子信号的到达太阳系质心SSB处的时间，令n的初始值为0；

步骤S2：根据脉冲星相位预测模型计算得到t_n时刻对应的脉冲星相位φ(t_n)；

步骤S3：由步骤S2得到的t_n时刻对应的脉冲星相位φ(t_n)，采用反函数法递推得到第n+1个脉冲星光子信号到达SSB处时对应的脉冲星相位φ(t_n+1)；

步骤S4：将步骤S3中得到的第n+1脉冲星光子信号对应的脉冲星相位φ(t_n+1)，代入脉冲星相位预测模型得到方程，求解该方程得到第n+1个脉冲星光子信号到达SSB处的时间t_n+1；

步骤S5：n的值加1，判断t_n是否小于模拟结束时刻t_f，当t_n<t_f时，转步骤S2；当t_n≥t_f时，结束本次模拟过程，此时得到脉冲星光子信号时间序列{t₀,t₁,t₂,…t_n,t_n+1,….t_f}，即实现了X射线脉冲星信号模拟。

2.根据权利要求1所述的一种高精度的X射线脉冲星信号模拟方法，其特征在于在：脉冲星相位预测模型为

$\mathrm{\&phi;}\left(t_{n+1}\right)=\mathrm{\&phi;}\left(t_{\mathrm{ref}}\right)+f[t_{n+1}-t_{\mathrm{ref}}]+\underset{m=2}{\overset{m=N}{\mathrm{\&Sigma;}}}\frac{\frac{d^{m-1}f}{{\mathrm{dt}}^{m-1}}}{m!}{(t_{n+1}-t_{\mathrm{ref}})}^m$

其中，N是脉冲星相位预测模型的阶次，N取值范围为[2，+∞)，f和分别是脉冲星的自转频率和自转频率的m-1阶倒数值，φ(t)是t时刻脉冲到达SSB处的相位，φ(t_ref)是t_ref参考时刻脉冲到达SSB处的相位，t、t_ref和t_n都采用太阳系质心力学时TDB作为时间系统。

3.根据权利要求1所述的一种高精度的X射线脉冲星信号模拟方法，其特征在于在：反函数递推公式如下：

φ(t_n+1)＝Λ^-1(Λ(φ(t_n))-ln U)

其中，U为[0,1]区间上均匀分布随机变量，Λ(·)为X射线脉冲星的累积泊松速率函数，Λ^-1(·)为累积泊松速率函数的反函数。