技术领域
[0001] 本
发明属于X射线脉冲星自主导航技术领域，尤其涉及一种基于贝叶斯估计的X射线脉冲星导航TOA估计方法.
背景技术
[0002] 目前，在天文导航中，X射线脉冲星导航是一种依赖脉冲星
信号源的新导航技术。在众多天体中，X射线脉冲星具有良好的周期
稳定性、流强大、旋转
频率可预测等特点。正因为脉冲星有上述优点，才被纳入到目前的
导航系统中作为深空探测任务的“灯塔”。TOA估计是导航系统的基本观测量，从而获得
航天器的
位置、
姿态等重要信息。
[0003] 国内外的高校和研究单位针对脉冲星达到时间进行了大量的研究工作。X射线脉冲星TOA估计研究的方法大致分为两类：一类是时域方法，一类是频域的方法。时域法的
精度受制于轮廓取样频率，频域的方法不受到
采样频率的限制。时域的方法更多样，而频域的方法非常有限。在国外，最大化似然函数、利用最小二乘估计、互相关技术等一系列方法求解TOA。利用
栅格化搜索的数值
算法完成似然函数的最大化。为了进一步提高精度；将脉冲TOA估计问题重构为一个循环移位参数估计问题，利用离散傅里叶变换完成似然函数的最大化。
[0004] 目前的算法在满足导航系统要求方面仍存在
缺陷，效率和精度都无法满足脉冲星导航系统的要求。利用累积轮廓模型求解TOA的算法受制于轮廓取样频率，无法获得高精度。利用
光子流量模型求解的算法较少，而且易受观测时间的限制，在较短时间内可以获得较高精度，否则无法保证系统时延精度的要求。
[0005] 综上所述，
现有技术存在的问题是：目前的算法在满足导航系统要求方面仍存在缺陷，效率和精度都无法满足脉冲星导航系统的要求。利用累积轮廓模型求解TOA的算法受制于轮廓取样频率，无法获得高精度。利用光子流量模型求解的算法较少，而且易受观测时间的限制，在较短时间内可以获得较高精度，否则无法保证系统时延精度的要求。
[0006] 解决上述技术问题的难度：
[0007] 时域法受制于轮廓取样频率，无法获得高精度。直接利用光子序列求解，计算量大，耗时长。因此，目前在有效观测时间内获得较高精度的TOA十分困难。
[0008] 解决上述技术问题的意义：
[0009] 国内外都已开展X射线脉冲星自主导航研究，而TOA精度直接影响着自主导航系统的精度。因此，提高TOA精度相当于提高导航的精度，可以为深空探测提供更精确的导航数据，推动航天事业的发展。
发明内容
[0010] 针对现有技术存在的问题，本发明提供了一种基于贝叶斯估计的X射线脉冲星导航TOA估计方法。
[0011] 本发明是这样实现的，一种基于贝叶斯估计的X射线脉冲星导航TOA估计方法，所述方法包括：
[0012] 第一步，光子计数率统计情况的整体趋势符合泊松分布的条件下，X射线光子到达时间序列可建模成非齐次泊松过程；PSR B0531+21脉冲星的流量特性符合泊松分布，建立泊松分布信号模型，分为时频模型和稳频模型两个部分；
[0013] 第二步，光子序列的信号模型进行傅里叶变换后转换到频域，得到带有时延估计参数的光子流量概率函数表达式；
[0014] 第三步，利用贝叶斯定理将流量概率表达式转
化成能够计算时延参数的似然函数求解；
[0015] 第四步，利用贝叶斯估计计算工具多模式嵌套抽样算法，进行
迭代，进而计算似然函数的参数估计值。
[0016] 进一步，所述基于贝叶斯估计的X射线脉冲星导航TOA估计方法的脉冲星泊松分布信号模型建立包括：
[0017] (1)选择观测时间间隔(t0，tf)，则Tf＝tf-t0，定义ti为第i个光子达到的时间，任意递增到达的时间集为{t0，t1，...，tp，tf}；
[0018] t0≤t1≤…≤tp≤tf；
[0019] 其中，tp为集合中的随机到达时间；
[0020] (2)设在{0,t}时间间隔内探测到的光子数目为Nt。λ(t)为单位时间内探测到的光子的数量；
[0021] 假设t0＝0，N0＝0；
[0022] 对于任意的0≤t1≤…≤tp≤tf，有N(t1)，N(t2)-N(t1)，...，N(tn)-N(tn-1)相互独立；
[0023] 在一个时间间隔ht内探测到一个光子事件的概率为；
[0024]
[0025] 在时间间隔ht内探测到多个光子事件的概率为；
[0026]
[0027] 那么这个随机事件{N(t)，t≥0}就是一个非齐次泊松过程，Nt为随机过程的变量；在时间t内，将在(0，t)时间内探测到的光子的数量为n的概率为P(Nt＝k)，则；
[0028]
[0029] Nt的均值和方差表示为；
[0030]
[0031] (3)随机选取的时间间隔(t,s)内可以探测到的光子数量记为Nt-Ns，Nt-Ns也是随机过程的随机变量，将在固定时间间隔内探测到的光子数量为n的概率为；
[0032]
[0033] 用 表示到达时间集，那么到达时间集的M维联合概率
密度函数为；
[0034]
[0035] 其中 表示泊松过程的联合速率；
[0036] (4)对于N点序列{x[n]}0≤N,它的离散傅里叶变换为；
[0037]
[0038] 在一个周期内探测到的光子数目为NT，则流量概率为；
[0039]
[0040] 其中X射线脉冲星的流量密度函数λ(t)≥0为总流量密度，由两部分组成，一部分自身的光子密度，另一部分是背景环境中的光子密度；
[0041] λ(t)＝λb+λsh(φdet(t))(ph/s)；
[0042] 其中h(φ)为周期脉冲轮廓；φdet(t)为探测
相位；λb表示的是来自环境中的有效背景流量，λs表示脉冲星射电源的流量；则探测相位的表达式为：
[0043]
[0044] 其中，fd(t)为多普勒频移，fs为X射线源频率，观测起点的初始相位为φ0∈[0,1]；
[0045] (6)在稳频模型条件下，探测器的速率v(t)＝v为已知常量，则稳频模型下的探测相位为；
[0046] φdet(t)＝φ0+f0(t-t0)；
[0047] 其中:
[0048] 流量密度函数：
[0049] λ(t；φ0,f0)＝λb+λsh[φ0+f0(t-t0)]；
[0050] 模型条件下，v(t)不是常量，fd(t)跟时间有关，则探测相位变为：
[0051]
[0052] 流量密度函数改写成；
[0053] λ(t；φ0,v(t))＝λb+λsh[φ0+fs(t-t0)+φd(t)]；
[0054] 在一个脉冲周期内，流量均值表示为；
[0055]
[0056] 光子到达时间的概率密度函数：
[0057]
[0058]
[0059] 进一步，所述基于贝叶斯估计的X射线脉冲星导航TOA估计方法的贝叶斯估计的算法求解TOA估计包括：
[0060] (1)计算一个周期的流量表达式。在一个周期内的流量表达式如下：
[0061]
[0062] 流量概率密度为；
[0063]
[0064] 改写成；
[0065]
[0066] 其中，h(φ0+f0T)为流量强度；
[0067] (2)计算在稳频模型下的流量概率表达式。在稳频模型M中，在不同观测起点，脉冲频率f0，观测起点的初始相位为φ0发生变化的情况下，贝叶斯定理得：
[0068]
[0069] 其中样本数量为k，初始相位φ0与观测到的脉冲频率f0为参数；
[0070] (3)计算贝叶斯定理计算流量概率的似然函数；在不同观测起点获取流量强度的值，再计算稳频模型下的流量概率。将流量概率函数代入到贝叶斯估计模型中，可得观测数据在傅里叶域的似然函数；观测数据在傅里叶域的似然函数如下：
[0071]
[0072] 其中，λ(t；φ0，f0)＝λb+λsh[φ0+f0t]；
[0073] (4)利用贝叶斯参数计算工具MCMC和MultiNest求解TOA估计；
[0074] 设置工具MultiNest的参数值，分别利用两种工具，在不同
信噪比的条件下，计算时延参数。
[0075] 本发明的另一目的在于提供一种应用所述基于贝叶斯估计的X射线脉冲星导航TOA估计方法的X射线脉冲星导航系统。
[0076] 综上所述，本发明的优点及积极效果为：
[0077] 在较短的观测时间内，脉冲星信号的信噪比较低时，基于贝叶斯估计的TOA算法的精度较高，说明该算法符合脉冲星导航系统的实时性高的要求。
附图说明
[0078] 图1是本发明
实施例提供的基于贝叶斯估计的X射线脉冲星导航TOA估计方法
流程图。
[0079] 图2是本发明实施例提供的脉冲星PSR B0531+21流量时变特性示意图。
[0080] 图3是本发明实施例提供的脉冲星PSR B0531+21流量概率特性示意图。
[0081] 图4是本发明实施例提供的时延估计分布直方图(Nmcmc＝0,SNR＝-15dB)。
[0082] 图5是本发明实施例提供的时延估计分布直方图(Nmcmc>0,SNR＝-15dB)。
具体实施方式
[0083] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0084] 本发明提出了一种基于贝叶斯估计的TOA估计方法，从而在有效观测时间内提高TOA估计精度，满足脉冲星导航未来的工程发展要求。
[0085] 下面结合附图对本发明的应用原理作详细的描述。
[0086] 如图1所示，本发明实施例提供的基于贝叶斯估计的X射线脉冲星导航TOA估计方法包括以下步骤：
[0087] S101：光子计数率统计情况的整体趋势符合泊松分布的条件下，X射线光子到达时间序列可建模成非齐次泊松过程；PSR B0531+21脉冲星的流量特性符合泊松分布，则建立泊松分布信号模型，分为时频模型和稳频模型两个部分；
[0088] S102：光子序列的信号模型进行傅里叶变换后转换到频域，得到带有时延估计参数的光子流量概率函数表达式；
[0089] S103：利用贝叶斯定理将流量概率表达式转化成能够计算时延参数的似然函数求解；
[0090] S104：利用贝叶斯估计计算工具多模式嵌套抽样算法(MultiNest)，进行迭代，进而计算似然函数的参数估计值。
[0091] 本发明实施例提供的基于贝叶斯估计的X射线脉冲星导航TOA估计方法具体包括以下步骤：
[0092] 第一步，脉冲星泊松分布信号模型建立
[0093] (1)选择观测时间间隔(t0，tf)，则Tf＝tf-t0。定义ti为第i个光子达到的时间，任意递增到达的时间集为{t0，t1，...，tp，tf}，即；
[0094] t0≤t1≤…≤tp≤tf；
[0095] 其中，tp为集合中的随机到达时间。
[0096] (2)设在{0,t}时间间隔内探测到的光子数目为Nt。λ(t)为单位时间内探测到的光子的数量；
[0097] 假设t0＝0，N0＝0；
[0098] 对于任意的0≤t1≤…≤tp≤tf，有N(t1)，N(t2)-N(t1)，...，N(tn)-N(tn-1)相互独立；
[0099] 在一个时间间隔ht内探测到一个光子事件的概率为；
[0100]
[0101] 在时间间隔ht内探测到多个光子事件的概率为；
[0102]
[0103] 那么这个随机事件{N(t)，t≥0}就是一个非齐次泊松过程，Nt为随机过程的变量。在一定时间t内，将在(0，t)时间内探测到的光子的数量为n的概率为P(Nt＝k)，则；
[0104]
[0105] Nt的均值和方差可表示为；
[0106]
[0107] (3)随机选取的时间间隔(t,s)内可以探测到的光子数量记为Nt-Ns，Nt-Ns也是随机过程的随机变量，将在固定时间间隔内探测到的光子数量为n的概率为；
[0108]
[0109] 用 表示到达时间集，那么到达时间集的M维联合概率密度函数为；
[0110]
[0111] 其中 表示泊松过程的联合速率。
[0112] (4)对于N点序列{x[n]}0≤N,它的离散傅里叶变换为；
[0113]
[0114] 将(6)式代入到进行傅里叶变换，在一个周期内探测到的光子数目为NT，则流量概率为；
[0115]
[0116] 其中X射线脉冲星的流量密度函数λ(t)≥0为总流量密度，它由两部分组成，一部分自身的光子密度，另一部分是背景环境中的光子密度，即；
[0117] λ(t)＝λb+λsh(φdet(t))(ph/s) (9)
[0118] 其中h(φ)为周期脉冲轮廓；φdet(t)为探测相位；λb表示的是来自环境中的有效背景流量，λs表示脉冲星射电源的流量。则探测相位的表达式为：
[0119]
[0120] 其中，fd(t)为多普勒频移，fs为X射线源频率，观测起点的初始相位为φ0∈[0,1]。
[0121] (6)在稳频模型条件下，航天器以恒定不变的速度开始移动，脉冲频率f0并不随时间变化或者变化的幅度很小。探测器的速率v(t)＝v为已知常量，则稳频模型下的探测相位为；
[0122] φdet(t)＝φ0+f0(t-t0) (11)
[0123] 其中:
[0124] 将(11)带入(6)式中得流量密度函数：
[0125] λ(t；φ0,f0)＝λb+λsh[φ0+f0(t-t0)] (12)
[0126] 模型条件下，v(t)不是常量，fd(t)跟时间有关，则探测相位变为：
[0127]
[0128] 流量密度函数改写成；
[0129] λ(t；φ0,v(t))＝λb+λsh[φ0+fs(t-t0)+φd(t)] (14)
[0130] 在一个脉冲周期内，流量均值表示为；
[0131]
[0132] 则将(14)式代入到(10)中得光子到达时间的概率密度函数：
[0133]
[0134]
[0135] (7)时变频率模型条件下的概率模型较为复杂，本发明暂时考虑在稳频模型下的流量概率函数作为计算依据。因此，在后续参数计算过程中，采用式(17)求解流量概率似然函数的后验分布参数值。
[0136] 第二步，贝叶斯估计的算法求解TOA估计
[0137] (1)计算一个周期的流量表达式。在一个周期内的流量表达式如下：
[0138]
[0139] 流量概率密度为；
[0140]
[0141] 形式可以改写成；
[0142]
[0143] 其中，h(φ0+f0T)为流量强度。
[0144] (2)计算在稳频模型下的流量概率表达式。在稳频模型M中，在不同观测起点，脉冲频率f0，观测起点的初始相位为φ0发生变化的情况下，式(20)代入贝叶斯定理得：
[0145]
[0146] 其中样本数量为k，初始相位φ0与观测到的脉冲频率f0为参数。
[0147] (3)计算贝叶斯定理计算流量概率的似然函数。在不同观测起点获取流量强度的值，再计算稳频模型下的流量概率。将流量概率函数代入到贝叶斯估计模型中，可得观测数据在傅里叶域的似然函数。在计算边缘后验分布时，不存在对贝叶斯估计的依赖，因为它跟参数本身的数值无关，则观测数据在傅里叶域的似然函数如下：
[0148]
[0149] 其中，λ(t：φ0，f0)＝λb+λsh[φ0+f0t]。
[0150] (4)利用贝叶斯参数计算工具MCMC和MultiNest求解TOA估计。设置工具MultiNest的参数值，分别利用两种工具，在不同信噪比的条件下，计算式(22)中的时延参数。
[0151] 下面结合具体实施例对本发明的应用原理作进一步的描述。
[0152] 本发明构建X射线脉冲星信号模型的过程：我国发射的X射线脉冲星试验卫星01星实现了对PSRB0531+21脉冲星的在轨观测，并获得该脉冲星的
辐射轮廓以及大量的科学观测数据。XPNAV-1卫星的主
载荷是北京控制工程研究所研制的掠入射聚焦型X射线脉冲星望远镜(Grazing Incidence Focusing X-Ray Pulsar Telescope,i FXPT)。根据MJD57727-57741期间的观测数据分析脉冲星信号的特性。
[0153] 根据观测特性绘制在观测过程中脉冲星PSR B0531+21光子计数率随时间变化的流量图。每个时间段的脉冲星PSR B0531+21光子计数率在平均流量值附近轻微浮动，图2中光子的平均流量值为14.7ph/s。
[0154] 对PSR B0531+21脉冲星的光子计数率的总体情况进行概率统计，与泊松分布的曲线进行拟合。如图3所示，光子计数率统计情况的整体趋势符合泊松分布。建立泊松分布的信号模型：样本数量为k，初始相位φ0与观测到的脉冲频率f0为参数。X射线脉冲星的流量密度函数λ(t)≥0为总流量密度，它由两部分组成，一部分自身的光子密度，另一部分是背景环境中的光子密度，流量概率的计算模型为；
[0155]
[0156] 获取脉冲星数据：采用RXTE航天器对脉冲星Crab的观测数据进行计算。可从美国高能
数据中心FTP上下载，下载路径为ftp://legacy.gsfc.nasa.gov/rxte/data/archive/。PSR B0531+21的周期为33.5ms，X射线光子流量为1.54ph/s/cm2，背景光子流量为0.005ph/s/cm2，光子序列的个数N为225，算法采样点数为112，代入到步骤4的公式中，得到含有时延参数的流量似然函数的模型。
[0157] 使用Matlab R2017a
软件读取实测数据，利用工具MultiNest计算TOA估计值。其中“Nmcmc”参数值表示抽样方案。本发明设置了两种抽样方案。如果设置了可选的“Nmcmc”参数并且Nmcmc>0，则将使用MCMC从方案中提取新的样本。这两种方法一旦达到容限标准，采样就会停止。本发明分别采用MultiNest和MCMC两种参数抽样方法得到TOA估计值。
[0158] 在Nmcmc＝0时，用MultiNest工具计算TOA估计，在信噪比为-15dB条件下，所求的TOA估计分布直方图如图4所示。经过1747次迭代运算后，最终得到的估计精度为1.68*10-5S。以此类推，设置不同信噪比，经计算得到的脉冲TOA的实验结果如表1所示。
[0159] 表1
[0160]信噪比(dB) -15 -10 -6 4 10 16 20
估计精度(×10-2s) 0.00168 0.00167 0.00159 0.00130 0.00129 0.00129 0.00128迭代次数 1747 1461 1373 928 886 763 664
耗时(s) 3.151 2.826 2.712 2.615 2.560 2.382 1.619
[0161] 在Nmcmc>0时，使用MCMC算法的计算结果如表2所示。同样地，信噪比为-15dB时，TOA估计分布直方图如图5所示。
[0162] 表2
[0163]信噪比(dB) -15 -10 -6 4 10 16 20
估计精度(×10-2s) 0.00200 0.00189 0.00179 0.00173 0.00169 0.00159 0.00132迭代次数 1630 1551 1345 1129 914 837 747
耗时(s) 4.926 4.653 4.167 3.736 2.56 2.398 2.498
[0164] 为了进一步校验算法性能，将本发明的算法与其他两种常见时域时延估计算法比较。分别在两种不同的观测时间长度下对比估计精度，实验结果如表3所示。
[0165] 表3
[0166]
[0167] 以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何
修改、等同替换和改进等，均应包含在本发明的保护范围之内。
